{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 3: Introduction to Probability Theory (Part I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "\n",
    "+ To understand the interprtation of probability theory as a representation of our state of knowledge.\n",
    "+ To learn about the common sense assumptions that give rise to the basic probability rules.\n",
    "+ To gain experience with the basic probability rules."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Readings\n",
    "\n",
    "+ These notes.\n",
    "+ [How random is a coin toss, Persi Diaconis, Video by Numberphile](https://youtu.be/AYnJv68T3MM).\n",
    "+ [Should you catch a tossed coin?, Persi Diaconis, Video by Numberphile](https://youtu.be/Obg7JPd6cmw).\n",
    "+ [Chapters 1 and 2, Probability Theory: The Logic of Science, E. T. Jaynes](https://bayes.wustl.edu/etj/prob/book.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('talk')\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability as a representation of our state of knowledge\n",
    "\n",
    "Let's call $I$ all the *information* you have a this given moment.\n",
    "And I am talking about absolutely everything: what your parents taught you, what you learned in school, what you learned in college, what your eyes see right now on some scientific instruments.\n",
    "Now consider a, well-defined, sentence $A$ that says something about the world.\n",
    "For example, $A$ could be \"The result of the next coin toss John performs will be heads.\"\n",
    "Or anything really.\n",
    "We want a technical machinery that can turn all the information $I$ we have into a real number that tells us how plausible it is that $A$ is true.\n",
    "This is what probability theory does.\n",
    "It gives us such a number.\n",
    "Call it $p(A|I)$ and read it as \"the probability that $A$ is true given that we know $I$.\"\n",
    "So, probability theory is an attempt to represent our state of knowledge about the world."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## But what about frequencies?\n",
    "\n",
    "In introductory courses to probability or statistics, we usually learn that the probability of an event is the frequency with each it occurs in nature.\n",
    "This is absolutely fine if the event is something that indeed occurs repeatedly.\n",
    "However, this intrepretation is quite restrictive.\n",
    "In particular, what can we say about an event that can happen only once?\n",
    "This interpretation forbids the quantification of epistemic uncertainties.\n",
    "We will expand the interpretation of probability.\n",
    "It can be shown, see E. T. Jaynes book for the proof, that this interpretation is compatible with the frequency interpretation.\n",
    "That is, when events occur repeatedly then the probabilities do become frequencies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Coin tossing\n",
    "\n",
    "Let's consider the classic coin tossing example.\n",
    "\n",
    "![Coin flipping](coin_flipping.png)\n",
    "\n",
    "The typical interpretation is:\n",
    "\n",
    "> The probability of heads = the observed frequency of getting heads in repeated random experiments\n",
    "\n",
    "But in what sense is tossing a coin a truly random process?\n",
    "As we will show below, it is not at all random.\n",
    "Newton's laws govern exactly what happens to the coin.\n",
    "The problem is that we do not know the initial conditions of the coin.\n",
    "This is a classic example of epistemic uncertainty misrepresented as aleatory uncertainty.\n",
    "Our intrepretation is:\n",
    "\n",
    "> The probability of heads = how sure I am about the outcome of a coin toss experiment given all the information I have at the moment\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the physics of an idealized 2D coin.\n",
    "Its dynamics can be described by the initial value problem:\n",
    "$$\n",
    "\\begin{array}{ccc}\n",
    "\\ddot{y} &=& -g,\\\\\n",
    "\\ddot{\\theta} &=& 0,\\\\\n",
    "y(0) &=& 0,\\\\\n",
    "\\theta(0) &=& 0,\\\\\n",
    "\\dot{y}(0) &=& v_0,\\\\\n",
    "\\dot{\\theta}(0) &=& \\omega_0.\n",
    "\\end{array}\n",
    "$$\n",
    "Here we have assumed that the coin is always toss vertically from position, $y(0)=0$, while laying flat, $\\theta(0)=0$, showing heads.\n",
    "$y(t)$ measures the vertical distance from the tossing finger, and $\\theta(t)$ measures the angle of the coin with respect to the horizontal axis.\n",
    "$v_0$ is the initial tossing velocity and $\\omega_0$ is tossing angular velocity.\n",
    "$g$ is the acceleration of gravity.\n",
    "\n",
    "Therefore, the only free parameters that determine the outcome of the coin are $v_0$ and $\\omega_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is trivial to solve this initial value problem. The solution is:\n",
    "$$\n",
    "\\begin{split}\n",
    "y(t) &=& -\\frac{1}{2}gt^2 + v_0t,\\\\\n",
    "\\theta(t) &=& \\omega_0t.\n",
    "\\end{split}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the first equation, we can find the time $t_1$ at which the coin returns to the hand. That time is:\n",
    "$$\n",
    "y(t_1) = 0 \\implies t_1 = \\frac{2v_0}{g}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this result, we can predict the angle the coin makes with the horizontal axis when it returns to the hand. That angle is:\n",
    "$$\n",
    "\\theta(t_1) = \\frac{2v_0\\omega_0}{g}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This angle is enough for predicting the outcome of a tossing experiment.\n",
    "Specifically, we will get **heads** if\n",
    "$$\n",
    "\\frac{2v_0\\omega_0}{g} (\\text{mod}\\;2\\pi) \\in \\left(0,\\frac{\\pi}{2}\\right)\\cup\\left(\\frac{3\\pi}{2},2\\pi\\right),\n",
    "$$\n",
    "and **tails** if\n",
    "$$\n",
    "\\frac{2v_0\\omega_0}{g} (\\text{mod}\\;2\\pi) \\in \\left(\\frac{\\pi}{2},\\frac{3\\pi}{2}\\right).\n",
    "$$\n",
    "Note that we have used $\\text{mod}\\;2\\pi$ to find a unique angle between $0$ and $2\\pi$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let $X\\in\\{H,T\\}$ be a variable describing the result of the coin toss. We have shown that, conditional on knowning $v_0$ and $\\omega_0$ we can predict $X$ exactly. Mathematically, we can write:\n",
    "$$\n",
    "X = \n",
    "\\begin{cases}\n",
    "T,&\\;\\text{if}\\;\\frac{2v_0\\omega_0}{g} (\\text{mod}\\;2\\pi) \\in \\left(\\frac{\\pi}{2},\\frac{3\\pi}{2}\\right),\\\\\n",
    "H,&\\;\\text{otherwise}.\n",
    "\\end{cases}\n",
    "$$\n",
    "Let's draw the corresponding graphical causal model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: coin_toss_g Pages: 1 -->\n",
       "<svg width=\"206pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 206.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>coin_toss_g</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 202,-112 202,4 -4,4\"/>\n",
       "<!-- omega0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>omega0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"18.44\" y=\"-86.8\" font-family=\"Times,serif\" font-size=\"14.00\">ω</text>\n",
       "<text text-anchor=\"start\" x=\"28.56\" y=\"-86.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">0</text>\n",
       "</g>\n",
       "<!-- X -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">X</text>\n",
       "</g>\n",
       "<!-- omega0&#45;&gt;X -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>omega0&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M41.57,-74.83C51.75,-64.94 65.52,-51.55 77.03,-40.36\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.47,-42.87 84.2,-33.38 74.59,-37.85 79.47,-42.87\"/>\n",
       "</g>\n",
       "<!-- v0 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>v0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"92\" y=\"-86.8\" font-family=\"Times,serif\" font-size=\"14.00\">v</text>\n",
       "<text text-anchor=\"start\" x=\"99\" y=\"-86.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">0</text>\n",
       "</g>\n",
       "<!-- v0&#45;&gt;X -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>v0&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-71.7C99,-63.98 99,-54.71 99,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-46.1 99,-36.1 95.5,-46.1 102.5,-46.1\"/>\n",
       "</g>\n",
       "<!-- g -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>g</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">g</text>\n",
       "</g>\n",
       "<!-- g&#45;&gt;X -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>g&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M156.43,-74.83C146.25,-64.94 132.48,-51.55 120.97,-40.36\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"123.41,-37.85 113.8,-33.38 118.53,-42.87 123.41,-37.85\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a22782d10>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Digraph\n",
    "gct = Digraph('coin_toss_g')\n",
    "gct.node('omega0', label='<&omega;<sub>0</sub>>')\n",
    "gct.node('v0', label='<v<sub>0</sub>>')\n",
    "gct.node('g', style='filled')\n",
    "gct.node('X')\n",
    "gct.edge('g', 'X')\n",
    "gct.edge('v0', 'X')\n",
    "gct.edge('omega0', 'X')\n",
    "gct.render('coin_toss_g', format='png')\n",
    "gct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation of a coin toss experiment\n",
    "\n",
    "A typical human throws the coin with\n",
    "$$\n",
    "v_0 \\approx 2.5\\;\\text{m/s},\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\omega_0 \\approx 200\\times 2\\pi\\;\\text{rad/s}.\n",
    "$$\n",
    "Let us simulate how sensitive is the result on the choice of these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.constants import g # acceleration of gravity\n",
    "# Let us code H as 0 and T as 1.\n",
    "def X(v0, omega0):\n",
    "    tmp = (2 * v0 * omega0 / g) % (2. * np.pi) # taking g = 9.8 m/s**2\n",
    "    if tmp > 0.5 * np.pi and tmp < 1.5 * np.pi:\n",
    "        return 'T'\n",
    "    return 'H'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H\n"
     ]
    }
   ],
   "source": [
    "# Try it out here:\n",
    "v0 = 2.5\n",
    "omega0 = 200 * 2. * np.pi\n",
    "print(X(v0, omega0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0d1159fb86be4d9e904f17c7faa540bd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=2.5, description='v0', max=3.0, min=2.0), FloatSlider(value=3453.31853…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's do some interactive tests\n",
    "from ipywidgets import interactive\n",
    "def print_X(v0, omega0):\n",
    "    print(\"X = \", X(v0, omega0))\n",
    "interactive(print_X, v0=(2., 3., 0.1), omega0=(100*2*np.pi, 1000*2*np.pi, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Play with the interactive tool above moving $v_0$ and $\\omega_0$ changing the result of the coin toss experiment. To which of the two parameters is the result of the coin toss most sensitive? \n",
    "\n",
    "+ Consider a coin tossing experiment in the moon. Would it be easier or harder to manipulate the result?\n",
    "\n",
    "+ To describe the dynamics of a real 3D coin you need 6 degrees of freedom. You need three degrees of freedom to describe the center of mass, two degrees of freedom to describe a unit vector pointing towards the instantaneous axis of rotation, and one degree of freedom to describe the angle of rotation. Using the Newton-Euler equations you can write down 6D dynamical system. How many arbitrary initial conditions do you have?\n",
    "\n",
    "+ List at least three other \"random\" experiments which are clearly deterministic. You can get inspired by gambling applications.\n",
    "\n",
    "+ Consider the example of a casino roulette. Discuss with your neighbors and report here. How many degrees of freedom do you need to describe its dynamics? What are the initial conditions that you need to know and what can you use to measure them?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rise of uncertainty\n",
    "\n",
    "As we saw, the outcome of a coin toss experiment is:\n",
    "$$\n",
    "X = \n",
    "\\begin{cases}\n",
    "T,&\\;\\text{if}\\;\\frac{2v_0\\omega_0}{g} (\\text{mod}\\;2\\pi) \\in \\left(\\frac{\\pi}{2},\\frac{3\\pi}{2}\\right),\\\\\n",
    "H,&\\;\\text{otherwise}.\n",
    "\\end{cases}\n",
    "$$\n",
    "Let's think of it as a function of $(v_0,\\omega_0)$.\n",
    "It is obvious that $X$ flips whenever\n",
    "$$\n",
    "\\frac{2v_0\\omega_0}{g} = 2\\pi k + \\frac{\\pi}{2},\n",
    "$$\n",
    "for $k=1,2,\\dots$, as well as whenever\n",
    "$$\n",
    "\\frac{2v_0\\omega_0}{g} = 2\\pi k + \\frac{3\\pi}{2}.\n",
    "$$\n",
    "Let's plot some of these curves in the $v_0-\\omega_0$ plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\omega_0$ rad/s')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEdCAYAAABQXlN8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd5gUxdaH39olo+wqIJj1iihcUTGBogICJsxZzOneK6AgSjJnRFRUUIJkJYqgSM4555zTAruwOaeZOd8fp6Z32A9UdoFdoN7n6Wdmuqurq1vs355Qp4yI4HA4HA5HcSesqAfgcDgcDsc/wQmWw+FwOE4InGA5HA6H44TACZbD4XA4TgicYDkcDofjhMAJlsPhcDhOCIqNYBljzjPGfGeMmWuMSTPGiDGmQb42ZxtjPjPGLDTGxBtjko0xS4wxzxljwvK1fd72caitzCGu38wYs8oYk2WM2WOM+eJQ7RwOh8NRNJQo6gGEUA14ElgOTAPuO0Sba4FngEHAR4AfuBcYANQGWh/inGeBLfn2ZYf+MMY8DfwM9LB91AA6AxcBTxTgXjDG+NA/CFIKcr7D4XCcglQAAiJySG0yxWXisDEmTEQC9vsDwGigoYjMDGlzBpAmIrn5zu2PClklEUmy+54H+gO1RWTlX1w3HNgDLBaR+0P2vwL0BuqKyKIC3E8AMBEREUd6qsPhcJySJCcnA4iIHNL7V2wsrKBY/U2bxMMcWgI8D5wNJB3hpesCVYGB+fYPBn4AHgaOWLCAlIiIiIikpCMdjsPhcJyaREZGkpycfFivVLGJYRWS24B0YOchjk00xvhtzGuYMaZavuNX2M+1oTtFJAPYFnLc4XA4HEVIsbGwCoox5kHUCvpYRDJDDsUAnwELgVTgGqAjsMgYc72IbLftKtrPhEN0nxByPP91/850cr5Ah8PhOIqc0IJljKmLJktMBT4JPSYiE4GJIbtmGWOmAMuAt4GX83V3uGBe8QjyORwOxynOCStYxpjrUUFaAdwvIr6/O0dE1hpjlgE3huyOt58VQ74HORPYcZi+Iv9mfEk4K8vhcDiOGidkDMsYcy0wGdgA3G3jTf+UMCA0wWOd/TwoVmWMKQdcQr7YlsPhcDiKhhNOsIwxtYEpaELEnSKSegTnXoHO11oYsnshGu96Jl/zJ4GSwKhCDdjhcDgcR4Vi5RI0xjxiv15vP+sbYyoB6SIywRhzGSpWAeB9oIYxJrSL9SKSYvuaAkxHLag0VKjaA8loMgYAIuIzxnQABhhjugMjyZs4PFJEQsXN4XA4HEVEsRIs4Nd8vz+0n7vQqhM3kpe1N+4Q5zcEZtrva4GngfOBskA0MAbNJtwdepKIDDTG+FFBewWIA3oCHxT4TgrDvHlw7bVQxlWGcjgcjiDFptLFyYYxJqnAE4fnz4fERGja9OgPzOFwOIopduJw8uGS2k64GNYpQd26sGBBUY/C4XA4ihVOsIojYWFgDPj+NlPf4XA4ThmcYBVXbr0V5swp6lE4HA5HscEJVjFk8WLIuakBzJhR1ENxOByOYoMTrGJIVhZMm10S/H4I/G0Re4fD4TglcIJVDKlXD+bOBW64AZYsKerhOBwOR7HACVYxJDxc8y5yG94OkyYV9XAcDoejWOAEq5jSsCHMXFQWMjPBzZVzOBwOJ1jFklWruPUmH7NmAVdeCatWFfWIHA6Ho8hxglUciY+nxBzNEPTdfR+MGVPEA3I4HI6ixwlWceTWW2HWLJ2Ktbw8ZGS4bEGHw3HK4wSrOFKiBBhDw5tzmT4duOUWmzbocDgcpy5OsIorjRtTcvY0AgHwN7odJk8u6hE5HA5HkeIEqxiyYQP4b7wZ5sxR7+D8klpbMCenqIfmcDgcRYYTrGLInj0wc044lChBo1tymDYNuOMONyfL4XCc0jjBKoY0bGjLCDZpQokZUyhZEjJr36TrZDkcDscpihOsYkiJElrpIuc6Fan77oMxY8OgXDlISyvq4TkcDkeR4ASrmHL77TBpShiULEntGlmsWAE8/DCMHFnUQ3M4HI4iwQlWcWT7dm6qG1APYNOmmHFjqVQJYivX1IwMh8PhOAVxglUc2biRsHlzKFsWMmpeB0uW8NhjMGIEcOmlsGlTUY/Q4XA4jjtOsIojTZrAlCnccw+MHWfgrLO4oPR+du+GPOVyOByOUwsnWMWRkiUhLIzaNbM1dvX44zB8ONWrw6boClrBPTe3qEfpcDgcxxUnWMWVu+7CTJxARAQknXYe7N3LY4/B8OHA3XfD+PFFPUKHw+E4rjjBKobs3QtSpy4sXMgjj1gP4JVXcvr2VWRnQ8719WDevKIepsPhcBxXnGAVQ5Ytg/kLDFSoQLXKyWzbBjz4IIwezQMPwB9jDFStCtHRRT1Uh8PhOG44wSqG3HknTJyIzrv67Tdq1oR1O8pBTg7XXZXLkiXA00/DL78U9VAdDofjuOEEqxhSqpRWusi68DLYvJlHHrHzhe++GzN+HBdfDNvTzoKEBJd84XA4ThmcYBVT7r0Xxo4FzjmH8ol7yMmB3Bs0dtWsGQwdCuof/KOoh+pwOBzHhWIjWMaY84wx3xlj5hpj0owxYoxpcJi2zYwxq4wxWcaYPcaYL4wxZQ7RrooxZqAxJs4Yk26MmWOMuakwfR4X4uK49hpR19+TT8KQITRtCuPGG7j4YiLit5OeDrm1b4BFi4pkiA6Hw3G8KTaCBVQDngTSgGmHa2SMeRoYDMwD7gI+B1oAA/K1K2P7qQ+8BjwIpALTjDG1C9LncWP2bMySxVSuDAekMiQkcGOdAAsWAE89BUOGcO+9MOZPA1dcAWvWFMkwHQ6H43hSnARrtoicJSJ3Av0P1cAYEw50AcaISHMRmSEiPwJtgMeNMXVCmr8I/Bt4SESGishkVLSiUUEqSJ/Hh3vugbFj84paNGqEmT6NKlUgJjMCMjOpe00OCxfiTSp2OByOk51iI1giEvgHzeoCVYGB+fYPBnKBh0P2PQisEZHlIdfIBoYCTYwxpxegz+NDqVIQHs4FlTOJigIaNYJp03jqKfj5Z+ChhzCjR3HVVbByYxkoUwaSko77MB0Oh+N4UmwE6x9yhf1cG7pTRDKAbSHHg20PamdZDYQDNQrQ5/HjgQfg99/5979h7fowOPNMqoTFEh8PvquuhWXLePRRa4E9+ywMzK+3DofDcXJxoglWRfuZcIhjCSHHg20P1y60ryPp08MYk/RXGxDxVzfyV6SlAVdfDStX5lW6aNYMhgzRxRzHADVqUHrbeiIi4ECZC3QScU5OQS/pcDgcxZ4TTbCCyD/cf7h2R9L2r/o4JowbB8uXA+efT7n4KIyBjDPPg337uLGuaPLFE0/A0KE8/7w1rp54wsWyHA7HSc2JJljx9vNQVs+ZHGwlxf9FO0LaHkmfHiIS+VcbkPxXN/JXPPAAjB6NitCwYXl5FfXqYebN5fLLYcOuclCuHFXC40hMhJyaV8OqVSDHXV8dDofjuHCiCdY6+3lQXMkYUw64hIPjUOvyt7PUAvzAxgL0eVwoXRrCwyGjXCWIi6NmDWH9eqBpUxg7lieftBOHX3gBBgzIcxvecQdMmnS8h+twOBzHhRNNsBYCMcAz+fY/CZQERoXsGw3UMsZcHdxhjCll204VkZQC9HncsGUE4dZbYdYsDWmtCYdzzqFcfBTh4ZBaviokJHDNv3XdLGnUGKZOLYrhOhwOxzGnWAmWMeYRY8wjwI12V3277y4AEfEBHYAHjDHdjTENjDGvAl2BkSKyMKS7vsAGYJQx5nFjTBNUfM4B3gk2OsI+jw85OdSqBWvXopVwJ0zIE7Bnn4VBg3j6aZvibv2FXiWMa66BpUuP+5AdDofjWFOsBAv41W7t7O8P7e8ewQYiMhC1huoDk4B3gZ7As6EdiUgWcBtavaIH8AcQCTQRkWX52v6jPo8bv/4Kq1ZRrRps3hYOVapQJmEfJUpAaokzICuLS87JZOdO8P37Kli9moYNhOnTIW+VR4fD4Ti5MOKC9McEY0xSRERERFJBJvRmZsLnn5PS9hO+/x7e/W8s9OvH1ofbM2UKvFp/PSxYwKIrXmLnTnj89PFQujSjkhtRuTLcsrW/WlpXXXXU78vhcDiOFZGRkSQnJyfbxLX/R3GzsBwAZctCeDgVwtPJzITcyMqQlES1C3PZvh0Cl9eEDRuoc4OwaBHIHXfC+PHBucZurSyHw3FS4gSruGJjU/fdZ1cQuf9+GDOGe+7ReVrUrw8zZ9KgAcyaEwZ16hC2aAG1asGq9SWhRg1YvbqIb8LhcDiOHk6wiiGBACo4GzdSp45dQaROHViwIJg06KW42zq5Xlphs2bWuHrmGZuV4XA4HCcHTrCKIUOGoPOurroKVq7k6qthxUoDNWpg1q/T3avD4MorCVu9kho1YN3GcKhVi1LrVlCzJqxcVxIuv9wtPeJwOE4anGAVQx56CIYNw7OaHn0URo5E6wkOHpxX+cLOIPZS3O1SxF4I65lnXFFch8Nx0uAEqxhSrhyEhUFqbhkoUYJSOWmcfjrEpZeFyEhKxUdTsSJEx5eCCy6gdNRWLrwQNm4rCdWqUXLLeq6+GpauLqXZggsWFPUtORwOR6FxglVMscXZCZpPzz0HAwYAL70E/frx4ovQpw/w/PMwYAAvvAD9+0OwoVe+KVgU101fcDgcJzhOsIojIlSvDlu2gPzrEti+nbOrBIiLg9wKFSE7mzNLpWEMxGeVh9NPp0xiNBdcAJt3lYbLLyd8zUquvx4WLAqDu+6C8eOL+q4cDoejUDjBKo4MGwabN3PzzTBvHgRTAR9+2MayrFX1yivw00/AK69A7968+CL064eWb/r5Z6/ohTS5HaZMAb+/SG/L4XA4CoMTrOLIfffBL7/kpazfeivMns3118OSJSAX/wt276ZKRR+ZmZBS4kwoXZqyifs4/3zYsqOEzstaMI+mTeHPsQaeegoGDy7qO3M4HI4C4wSrOFK+PISHUyIjhYoVIWa/8VYgvusumDgRePRRGD6cl1+2saz//Q969cqLbT3yCPz6K00aCzNmgK/29bBuHaSnF/HNORwOR8FwglVcsW69F16Avn3xito2bqzePa6/HpYt4/xz/MTHQ0apSDj9dMrGRXHZZbBiVRjcfjtMmpSXkPHaa9CtWxHfmMPhcBQMJ1jFlYsvhl27qHRmgKwsSMspBZUqYaL3ceONMH8+Wr5pxIg8q+q//4VevYIrkGiyxYQJXPlvPzt2QGrEeVCiBOzcWbT35nA4HAXACVYxZMQI2LYNXUF48mSef97O/7XK9NBDMGoUWq5pyRIuucjPgQOQIqdDxYqU2LOT+vVh2nTjJWi0aAHduwMtWsAPPxTl7TkcDkeBcIJVDLn7blu54rbbYNo0LrkEdu0C3+lnQHg44Ylx1Khhqy49+ij8+iuvvgo9ewL/+Q/07BmslYtcXRu2bePc01MoVQq27SurQjdjRhHfpcPhcBwZTrCKIaedBuHhkJRsoHp12LgxL6X9P/+B3r155hnr9rvxRli8mHPPVtdhXGZ5qF4ds2J5MO/Ci115xtXDD+s6JDk5RXujDofDcQQ4wSqmeG5AW3W9Th1YvBikUmXIzaVURpIWuV2JWlkjRtC8OfTogVa7GDiQW24WFi+GzMizoXRpysTspFEjGDfeQPPm8OOPRXuTDofDcQQ4wSqmnH8+7NsHvhJloFIliIqiUSOYNg2CM4afecYWub3xRliyhEoVcihdGvbGhHvVLTyrygaxmjaFqVMh68LLIDcXtm8v4jt1OByOf4YTrOLIqFGwe3ewWLvWD+zbl7vv1sUb5exzICWFEllp1K6tllewxuCrr1qBuuMOmDqVi8/34fNBVFxZuOUWmDiRli1tm9df1zR3V2fQ4XCcADjBKo40agT9+3PDDVrZggoVoGRJTHwcjRrB9OnAyy9Dr15ebVtq1oTduzmdVKpVg+UrjLoG+/fPm351330waRKXnJuFCGzZXVr3jRxZtPfrcDgc/4CjIljGmBLGmIeNMa8YY6oejT5PaSIidH2RxETq1YPZs1E3YJ8+wYWGkQsuhNRUwpMTqFsX5szBi0s995wmZMhVV8OuXZTPTuCaa2DOXKNW1fff5xlXDRrCihUQH1/Ud+1wOBx/yRELljHmS2PMkpDfBpgKjAB6AWuMMZccvSGeolgX3/33wx9/AGedBRkZmLTUPCvLZlkEMwjl3PPA5yM8Nob774fRo4E33oCuXXn8cc0YzD3vYggPp9TurTRrZhM72rWDzp2L9n4dDofjbyiIhXUnMCfk973ArUAXoJnd16GQ43Kccw7ExxOWk8X118OiRagbsHfvPCur8lkQHk5YzD4eeUQnHAeTKxo21Erv2adVhIsuwixfRvPmNnbVsiV060bdOsKOHRCTFQkNGmiqu8PhcBRTCiJY5wNbQn7fC+wQkQ4iMgzoCTQ6GoM7VRkzRjMEg2vd27nBcMEFkJKCSU6iUSNbU/DVV6FHD265RZMvsspEalmn5ct59VWbuW5z5C+vHiAnB7bvLa21CQcN4q234Kuv0NnKixY516DD4Si2FESwSgGhCys1RF2CQbYDZxdmUKc6t95qawPWrAkbNxKOn6uuguXL8QSqaVPNGAycHgGVK8PWrcGC7V45pmqXCNnZsGtPuIrfwIHBEBZyUz3YsYPTU/dxxx02caN9e/jiiyK9d4fD4TgcBRGsKKAugDHm38C/gFkhx88C0go/tFOXyEitdBEfDzz0EIwcmbfkfdWq4Pdj4mJp1gyGDMErx3TppZCYCHGJ4VoYd8gQWrWC774DbrgBduygTPJ+HnjAug/btoUuXWjSWFi3DvZlRGqF96FDi/T+HQ6H41AURLCGAc8ZY8YCY4EUIHT99drAtqMwtlOal1+2qwnfdBMsWkQJ46dGDVi1CoKFA+vU0d9ZlIG6dWHGDF57zQpUvXqwejVlfancfrtN3HjrLejShQYNYPVqOJBeXlczHj6c9u0170IaN4GtW2HHjiK9f4fD4chPQQSrEzAAuBEQ4FkRSQIwxkQA9wHTjtYAT1WqVIGMDEhNRa2l4cN55hmb1VexIpQqBfv28d//Qu/eePUBK0b6uegiWLYMaNUKvv2WO+/UtPe0sArQsCGMGZMnULc1gnXrKJ+8jyeegH790KzBb74Bn68oH4HD4XAcxBELlohki8hLIlJRRP4lImNCDqei8asPj9YA82OMGWCMkb/Yqtp2Mw9zfNgh+jzNGPO9MSbaGJNpjFlqjLnvWN3D3zJuHOzfHyxwodXVly2jpPFx4412XpYtaVGtmroOExINwZUaX3hBhS1Q9RyNb61YwRtvqAbRtCnMnUsFSeauu6xrsGNH6NyZG+sEiI2FjTtKa8p8165F9ggcDocjP/9IsIwxUcaYbsaYRsaY8MO1E5GAiCSLSO7RG+L/4xPUugvdbgUygUUiEhPSdssh2r57iD5HA0/ZY02B9cBoY8zdx+ge/pq6daF3by68EPbvh8xMCAaxHnlEyzXJ6RU0G3DlSlq21EQKrr4adu4kLCUpWORCJxz368e5Vf2cey4sXIjGrr78ksaNYe1aiEkpp0V2e/Xirbd0QnHWxTXUkps5s0gegcPhcOTnn1pYY4AHgClArDHmZ2PMg8aYcsduaIdGRLaJyMLQDagMlAX65muekb+tiGwNbWBFqTHwsoj0FZHpwHPAAuDr43BL/5+KFbXSxYEDvPSSzRi87jpYtQrj93HffTYmZbMBK1cSqla162O9/jp8+y3XXqtraMUlhmtSRo8evPiiJmlknV4ZrrkGxo/3EgPl2usgM5MSm9bRrh106oRabOPHQ0zMXw7X4XA4jgf/SLBEpIWInI9mB/YCrgN+A+KMMX8YY543xlQ8huP8O14EMoDhBTj3QSAZ+CO4Q0QEGAhcboypeVRGeKTY4FS1amplpaSgVtCAATRqpIaPjxLq4vvzT15+WYVNKp+lpd6XLuWNN+Drr4FatXT+1q6dtGljXYMPPwwzZ3JaVhyPP24TPFq1gp49ufCsTK66Cv4YY+CDD+Czz1w8y+FwFDlHFMMSkcUi0lFEagA1UfdcVdSyiTHGzDDGvG6MufAYjPWQGGPORqtvjBSRlHyHLzPGJBpjfMaYLcaYd40xJfO1uQJYLyKBfPtXhxw//lSqpJ9xcXmrCV91FWzbBmlpPPccDBgANGkCM2dSIpDDI4/AsGGoZTRoEGec7uO662DyZAgGsS66UKhc2RbVfftt6NSJG+sKqamwdkO4zsX67DMeegiWLoWdseU1ntWlS5E8BofD4QhS4OK3IrJRRDqJSB3gAuANdELxV8B2Y8xyY8ydR2mcf8VzQDj/3x04x47pQdSdOQv4GK15GEpFIOEQ/SaEHP9/GGOS/moDIgp2Oyow8fGoldWrF+eei64mHAfBBa5q19bM88RE1OXXqxe33KILOqamh2nsqlcvHn4YJk2CVF9ZaNYM+vblpZfg558ho1SkWmhDhtC6tU46zqx4nqbEDxvGu+9qFYysi2tAtWrWD+lwOBxFw1Gp1i4ie0Wku4g0RicOvwDs5PhYJ88DW0Vkdr4xvSciPUVkpoiMFZGX0ZT8B4wxN+fr468WhDrui0XVrm1XDq5cGQIBiIvLWyD4vPN03969tG5t51xdfrn6DKOiaNXKJvfVqgXJyRAVFZwfrMkcBw4QtmMb7dvbONVtt8GOHYTv3Ea7dvD55+jij9u3U3rnJjp0gI8/Rlc1XrsW1q073o/D4XA4gGOwHpaIJInIIBF5SES+Otr9h2KF5zKg/z88ZaD9vDFkXzyHtqLOtJ+Hsr4Qkci/2tC4WIGoXFlzLmJiIFixtlIlKFMG9uyB4AqMlStrEfd161CX33ffcc45Wjd3yRKgdWv45huqVhEuvdQuQfLmm9C1K+dW9XPNNVq3kHbtoGtXzj8rmzp1bJmmdu2gWzfOOyOdxo1txmGHDmqGJRzykTgcDscx5URfwPFF1A058O8aWoL3GxqvWgfUMMbkfxa17Ofagg+v4Hhxq4oVdQHH7dtp3twuxHj66fCvf8HKlbzyik4clrLloH59GD+eF19Ul19uyXKaXDFoEE8/renwqTml1YXYvTsPPqj1bqPjS2kVjE6duOceDZOt31wC3n0X3n+f2xoKaWmwYHG4mlsffOCSMBwOx3HnbwXLGDO9ANsxr3RhjCkPPApMEpG9//C0Z+3nwpB9o4FItOp8/rabRGR9oQZaQM44A047DaKiIFjV9rTT4Ior7FwqO0m4ZHiA++5TMeLee2HqVMKyMoIriMDNN8OuXZg9UcEcC7jySi1WuHSpt89//kWa1PHHH7Rvr+7H1PJVtar799/TsqWuPrI7JVItvA8/BDnu3lKHw3EqIyJ/uaGxqB35tgOolRJAXWaJIb8PANv/rt/Cbqh1JcDDhzh2CzDOtmmElovqa8c3Il9bA0wH4mz7hmjpqQBwbyHGlxQRESEFYto0kcRESUkRefddu2/gQJElS8TvF3ntNRG/X0QWLND9ItKunUhioojs2iXyySciIvLNNyIbNohIRobI66+LBALy558io0aJdvD66yLJybJqlciXX9rrfPihyJYtsn+/yBtviAQCItK/v8ikSZKVJdKihUhqqojMnCnSrVvB7s/hcDgOQUREhABJcrj36uEOHPYErc6+C+gKVA3ZXxX41grcxUfabwHGMQeIBUoe4lg1K1h7gCx0jtYKoDUQfoj2FYDuQIxtvxx4oJDjK7hgRUeLfPqpiIh8/73Ili0i4vOJtGolEgjIggUigwbZtm+/LZKQIAcOiLzzjt3Xu7fIokWSk6MCk5srIgsXivz0k4ioCO7dKyIHDoi0bSsSCMjAgSKTJ4t4J6WmyqJFKnoiIvLeeyKbNsn+/SKtW1vBHDJEZPTogt2jw+Fw5ONYCNYfwNC/OD4M+P1I+z3ZtkIJlohIp04ie/ZIVpbIm2/afZMmifzxh4ioRZWaKiKxsSomIvLzzyKzZ4t4ZlhmpqxfL9Kliz2/c2eR9eslPV0P+3yi1py10jp0EImKEpGYGL1oICBDhthL5uaKtGwpEhsrq1eLfPSR7fOrr0Tmzy/4fTocDofl7wSrIEkXDTh4/av8zLRuNUdhaN4cevSgdGm48UaYMQNdq2rWLMjM5PXXbfp6pUpw0UWwZAlPPaUrE2fnhhEsaVGjBpQvb7MG27SBH36gXFgWr7xiz7/tNjhwAFat4v33NZ6VHVlF09i7d+fJJ7Xk09qNJbTixQcfUKt6NjffbNPs27TRYr0u3d3hcBxjCiJYAtT4i+P/pgjmLp10VKgAZ58NGzfy0EOa8OD3A6+9Bt26ce65mv6+ciVaU3DQIExuDq1a2dJLF12kDRYv5r//tROFc0po4dtOnahVS5cwmTYNFZ3+/SmblUjbtjaf4oY6EBEB48fToYOWfTqQVUHT3d97j9saCmeeCb+NMvDRR1pWfteuInxgDofjpOdwptfhNrRSRA6aRWdC9hu06kQO+RIbTsWNQrgE58wRSUkRCfUHLl0q0q+fbdCli8i2beLzqZfO5xMNdHXuLCIivXppPobnGkxNlT17QmJcI0eKjBkjIiLt29t4VmKiJmH4/TJzpkjPnrbtp5+KrFolqakizZuLpKWJyJIlIp9/LiIi330nMn26HWuLFhoXczgcjgJwLGJY56GZgn5gH+oenGm/+9Gki/OOtN+TbSuMYO3ZoyEsEdHEhjlzRERjTMnJIpKZqSl8IrJqlciPP9q23bqJrFwpgYBqT3q6aAJHhw4iIvLnnyIjRti2774rsn27pKerzuTk2M6++EJEND9j2jRR0WvdWiQ6WqKjNe/D5xPN0LAX/vRTkUWLRDxVi48v0H07HI5Tm6MuWKIv4wjgc3TSbabd1tl9kQXp82TbCiNYImos7d4tmlduMySiozXrXERExo3zEjA++cS29fk89dmzx8vF0LZWqT7+2GYdZmeruGRkyNatIenzv/6qm+i+jRtFzapXXxVJS5P167XfQEBEhg0TGTZMAgG13tasEZGkJG2blFTge3c4HKcmx0Sw3HbsBSstTTMBRURk5UrPR9e7t8iyZXZ/27YiiYmSnu5Ns9KJV3ZS1bBhmlgoIj906xoAACAASURBVCLvvy+yc6eXtZ6ZKZoS2KGDSCAgU6aI9O1r2371lcjixZKbq1q5f7+opfbaayI5OTJnjjYRER3X+PHi94u89ZYVw/h4Fa2UlALfv8PhOPVwgnWCCpaIyIABIosX2x/vvScSG+uFpXw+0ZT2t98WEY1Z2WlWKiI21bxdO9UaycxUiyo7W3bt8k7TAJRVqu7dRebOFVW+t94S2b3b8/Klp4vI5s0a9AoEZPz4kDhX164i06ZJTo56D7dsEY1lNW9ufZgOh8Px9xwzwUIXcWyBLiv/fr7tvYL2e7JshRKshQtF0tO9YhSBgKjVYrMmDopbDR8uMmWKiGgexNatoie0aiWSnCypqSGJGVu3qqUlGoLykjh++EFk9mwJBFTgtm+XPIFLTpa9e1WIfD7RhAsbYBs+XGTwYNvHF1+IzJnjidbWrSISF6d9JCYW7Dk4HI5TimORdFEWmGATLAIhn6Hf/Ufa78m2FUqwdu70kh9mzRIZOtTu79fPZjdo3CoqSlSc3nhDJDVVsrNDxGn/fs+nuHp1yOThP/9UX6Fojsa8ebaPt98W2bZNsrPVZZiQIGrBtWghkp0t69eLdOxoxXPyZK8sU9++Ir/9Zvv4+GORBQskJ0f1cts2UbFq3lzFy+FwOP6CYyFYnawwfQzUtwL1DHCHzRZcBFx2pP2ebFuhXYJffmlNHQ0zJSZKXpp6drakp3uVmkT27TvI+ura1fYxcaJmGYoKy8yZdv/nn4usWyeBgIbB9uwR8ZQqKUkSEz2dEtmxQ1Pr/X5ZvFjks89sH2PGaEBNNI1+9GjRwXzyicjcuZKTI9Kmjcj69aJuwRYt7IUcDofj0BwLwdoCDLPfK1rBus3+LoHW7Ot0pP2ebFuhBSsjQ9/4ouEgL4tv82YvqWL+/JBEidGjvazBnj2t5SSiCrN+vQQCGn6KihI1wV57TSQhQTIyND8iPV3UCmrZUiQnR3bsUDHz+0WTPj76SCQQkGnTdO6ViGjmoS1q2KOHvXwgoII4a5b4fHrNFSskrwDv5s0FfyYOh+Ok5lgIVhbwqv0eaQXrzpDj7YAdR9rvybYdjaQLGT3aVqRVb+DChXb/t9+KrF0rIjoHyhpi6tbbt8/zEsbHS14NwNRUychQ71xmpqjJZsVp376QgrYhiRXLl3s6peaZTQ384w+1qkRE5JdfvEDWjz+GWFpffikyZYr4/Ro2mzdPdLJX27Yiy5cX7rk4HI6TkmMhWLHAa/Z7OOADXgw5/l8g80j7Pdm2wgjW8uVWVILJE9nZ4vervuTmSp4I+XySlaVf/X7RXHi79khiYsjyIDExmvUXCMju3Z4eqTjZyVqrV2sISkQ0TmYTK2bNCqnYPmGCJmiILlHSp4/dP3iwZ2n16RMSc+veXeS33yQQ0Bjan3+KDvSzz9Rd6XA4HCEci+K324DqACLiRycMPwJgjDHAQ0BUAfp1WCIi7OKLxhBcZjgsTFch7tYNKFHCWzW4dGl48UX44Qe0yu2zz8KPPxIZCY8/Dj16oEUDH3gAevbk/PPhrrt0pXsuvRQaNIDevalVC+rWhZ9+Am64QRdz/Oknbr0VqlXTUoHceSeccw7068eDD+o4+/cHmjWDUqVgwABeeklrHg4YALRoAWlpmIEDeOstrbE76JcwePtt2L0bBg0qoifscDhORAoiWFOBh40x4fZ3L+BOY8w2NL7VGF0s0VFA/vUvCAuDLVuA6tV1deANG6hZU/evWQPUqqX7ly+ndm39ungxcN11Kmjz5lGnDpQsCTNnAvXq6Y8pU6hfH0qXhgkT0GrtZcrAmDE0aQJly8KoUaiqlSsHw4Zx7736dfBgVPjKl4cBA3jkEV0VuU8fVB1POw169uSpp7R2b/fuqIBGRkLXrrz4ghAZCV99BfLyK7qs8hdfQCBQVI/a4XCcSBzO9DrcBpwGXAaUCNnXBl30cAnQnpCiuKfqRiFjWFlZIXOwguUpfD6v+lJ2tuRlDaanSyCgcajERNGT3nxTZP9+CQQ0y3DXLtuxXYhRRONfK1fa/SHrWnXv7k3t0tnIv/8uIrrelk061IQLm/Hxxx96joio2/DLL70EjY8+su7KOXM0xpabK4sWaSgrK0vU//nGG7aqrsPhOJU5qjEsdA7Ws0CdIznvVNwKK1giIjNmhMSDli/XVDzR0JMtzC6hRQPj4jSxMBAQLURrRS44BzgjQ/LiXwkJ4vdre28+19tv2+KBmuhnp3zpnCtb42nAAJ0wLCI6Acumtk+apHGqQEC0XMb774v4/bJihYpTdrZo36+/LpKWJjt3hkzPio7WH1FRhXpeDofjxOZoC1YYunzI/47kvFNxK5RgrV1ry6fr3F+vUMRnn3kpgT17hpRtGjXKS2mfMyckg2/dOjWjRHXN5l3orGBrpgUzBxMTRcWsVSuRXbskENBU+oMssKlTRUSkf/8QIf39d5HvvxcRXe34ww+tRbVqlaphZqZs3x4yGTkmxhOnpCQ1ENeuFc0yad/eLpnscDhORY5FluBWoN2RnneqbYUSrE2bvBTy2NiQIrjBNadyc70EQq8o+vvv29ISaojZFUl0gq9NO1+61JvCddCE4ISEkHqBmZlqgUVHi9+v7sQ1a+w5X3/tWVqDB6twiYgK2aefigQCsmqVdpuVJVpC3i43EtTIbdvsNd56S2ThQvH5VId//11UTXv21BsIBAr27BwOxwnLsRCs94A1QOkjPfdU2grtEvzmG1smQjXHhpHUHAkRszfftO/2YKWKjIyDJwmL6LwtG58aM0ZjUSKiM3ptLvu+fV4xdlUu66/z+7WvdetC+ho/XkTUsPNqGi5apLWbfD7ZuVOHkpws4pXN2L5dsrNVAL0Cu127ekGxYcPUzenFu9q2dXEth+MU41gIViO0msUG4DXgTuDW/NuR9nuybYUWrJBECxEVjdhYe6xnT28W8ezZIfOhoqJUNEQOniQcrGq7c6eIqAEzfbo9Z+pUL2Ni61YVQJ9PVG2sdeTz6emepdW9u6qVaI5F585WNNev18yP9HSJjdXTd+0S8ZRqzhwJBPR0zzobO1bkgw9EcnJk1SoVTS+u1bKlN0Ha4XCc/BwLwQrk2/z5Nlf89mgIlogmWlgTJikpJKEiuKSwXW/q6681ZCQimt5ny7Dv3q2GSiAg6qMLqZz+8cchBSdGjvQUZN06FSe/X7Ttq696ltbbb2uxdhHR7IsBA0REtdMmAKqp1ry5yIEDkpWl11+82I65e3dvgvHkyZorkpMjmkXSooVITIwkJanmLV4sqpydO+s5zkXocJz0HAvBeu6fbEfa78m2FUawtmzxDCsN8NjY1NSpWglJRA6qXuH3azzLW5n+hx+0RIXoi9+LWwUDSZmZntvQZrhrUMpmUqxZowaR3y95ltaBAxIIaJq6Fx8bPdpLuNi0SceQni6aofjaayKbNnlVLkaOtOeMH6+ZGTk5snWrDmffPnvOW2+JzJwpgYAmJnbvbnVq6lQ95lYxdjhOatwCjiegYK1ZE1JxPSTRQkQrJm3YYI/NneutohjMuMvNFX3Ld+xoF6XSBEIvbhUVpfOefD4vKXD3bnusb19bDFBLNXXsaEUruIrj3r2eAI0bZ8+ZNk0FyOfzjKt9+0QH8sEHXi3E337TVHmvXqG1qNLT1aKbMcOOu08fNRltdfhWrbT4r8TGqunlVfV1OBwnG06wTkDBEtHcBs/Nt2GDN/EqOI0qI8Me+/FHL6Fiw4aQeoDZ2aoe1irp0UO1RUQ0LvTuuyKBgJcUuG+fPdazpxefWr1aXXq5uaKm02uv2eWEVVcGDrTnrFmjYpKRIenpqofe2Pv10z4DAVm9Wj2ZSUmiIvjmmyLz50sgoFrZtasVtBUrvEElJ6ugTZwoKmg//aTKl51d4GfrcDiKJ06wTlDBOqiihYgqxNy5IqJGkrfEfSCgwa39+0XkoCx2zV4I6eTTT0NiUIsW6dpVgYCkp2uz6OiQa40YISJqDLVubdPUc3LUNbdihYioMfb119Ztt3evCuT+/eL367XGjLH9zZyp5lpmpsTFqe6tXGnH/uOP6v8LBGTFCj0WHS0qaB062Iq5mkX4/vtWqLds0QGvXl3g5+twOIofTrBOUMES0ffy55/bH8HaS9ZiGjcuRJhSUlQssrJERGM/M2bYYyHzrYKeQpstr8Eoa7kFRSsmxh4bONBLOd+1y1uhRE2gDz/0zLW5c1VXsrNFG7Ru7c02HjpU3Yd+v6jfsUULkZ07xe/XzHzPQgv6/uLjJS1N+5swwR4bOVKtwbQ02b1bBc1LyPjuO72AnWTtcDhObJxgnYiCtXu3fcurB8xLcoiLC1m4Sq0bz2LauVN9Z4GABAJqjXhzp1at0nhSICA+n7rsduywx6ZO9QJmaWmqKXv32mNDh3qll/bv9wwo5ccfvcDYtm16Xny86Ng6dfLcisuW6ZBTUkQF9Z13vKVFpkzRxI/UVNGMxDZtPKUdPlyzCNPTRct0vP66yPz54vfrkD77zOrzxo2qpm6NLYfjhOekEyygASCH2S7P17YJsBDIBA6gleUjD9HnacD3QLRtuxS4r5DjLLhgLVvmZd8FPX7eHKxly7wFqoILNXquvIULPfHx+dRo8Y7NnOkdy85WS8VLtpg8WU2eQMBbGNgmJmpWn51oFUz+s2Esnc1sjwXXg7SlCFVxvvhCxOeTAwf0PM+D98svqji5uRIdrceWL7c3NHCgHsvOlqgoPbZkieQlZHz6qVfu6bXXbPjO79cg3SefuMnGDscJzMksWO2Auvm2Mvna5QK/okuePGsFaR4Qlq/PKUA88BJwGzDIzim7uxDjLJxLsFs3raUkmln++uueYaUvdZt9F0zgs95AFYpffxWRvIIVCQn22NixXgHdYDV4r4r7zJnqfwwEJCdHvYjenN3589Uy8vkkJ0fdil6y3pIlqqjp6ZKbq0kftqyhKpSdCezzqffOTsPSDJEWLUSiosTvV+9e9+72HoMxqrVrPS3q1Mne444dOvB58yQQ0Olj779vLbioKB2LrcThcDhOLI6rYAHnAzWB849mv/muERSsB/6m3WJbkSMsZF8Te+7jIfvutvseDNlngLnAhkKMs3CCFVw6JDVVRNQC8VLdRTTrwvr1duwImSAsomJnY0yJiSpathvNlLClMf6faM2b51Wv9fn0Ep4wrV+vjVNTJRBQgfFiaMF8dtvR0KEhZZYOmgmsbsAOHex4MjJUbazALl2qVtquXaKpid98o+mSPp9s3arHPGtrwAB1c6akSHS0CuyoUfYZ/PGHCpd3Yw6H40TgmAsWUBL4yloowUoXAfu7C1CqsNfId72/FSzgXNumzSGO7QFGhPz+CUg6hNX1iu2jZgHHWfhKFzExGuSx9OkTUlIpOD/LJmEsWeKtaq988olnoe3f780XVkaM8CpbZGWp69Bz8y1frrGwnBwJBDROZsNReYGsPXtERLMAvblVwWrrdsLyypWqt7GxklflwmYDxsSo9nnxt/Hj9dzkZMnIUM38+WcrPqtXq1KtXy9+v1ai/+gj6/mLjlaltnPHJkxQbdy+XdS8/PxzvQHP/HQ4HMWZ4yFYvYH5wO1AJaCE/bzdut96FfYa+a4XFKz9gA9IBsYC14a0ucO2uf0Q508A1of8XgDMP0S7OraPxwo4zgIL1oEDIdbSxImeHy0QUM+cF18K5ojbScXjxnlrKkr+Uha7d6tIeO/uESO8hIqcHH3ve0uJbN2q/dp40C+/ePOT1Sp6801PDFet0n4TEiQvzvT99yKBgCQnq6Fjs/HVynr9dZHYWPH7tSBH1662qseBA9rYWoYzZ6qQ7tkjen9du2qczca2Wrf2cjf0i80kycxU665zZ5sCv3mzHhs92pV3cjiKOcdDsBKAioc5VhlILOw18vVZG+gK3A/cAvwH2A5kYReWBJpZsbnuEOcPBvaH/N4MjD1Eu0ttH68eZhxJf7NJQQVrwQJ9mXt8/bUnEMGauMnJ9timTapi9mXcv39IDCk3V1/WVuF27VId8iytP//0kjv8fvXOeS7AmBitI2jz3KdODan9FzS9bNp7QoLqkJdUsXixt/xxIKBWkSdMSUlqwdm89aABtXWr7XfIEL1Qaqqkpam11a+fvb3Nmw+KX40apdoZFSWqTp07ax2qjAwvm3/kSHvu5Ml6rqfKDoejuHE8BCseqHKYY2cDCYW9xj8YQ1UgDpgqBwvWtYdoOxiICfm9GfjzEO2CgnXIxSqPpWCJaNa4Jx7BOVgHDoiIutlatQqpNzhzpgaVLN99F7LEfU6OnmvjXcHMO69SxpQpms1n0+G/+srzsGmgqVUrr0z75s0qLl7Nwj//1CyL3Fzx+1UvvFqHwRT8BQtERC2xli1DljwZMULnV6WnS3a2eu969bLuxX37VGhtYsn8+TrmjRvtsxg0SANsdt7Wp59qqCs7W9Qf2KaN1oIKBGTaNL2FJUtEBfynnzSI5lY3djiKHcdDsL4DlgMPAhcDZ9rPB9H08G8Le41/OI4hQJr9fkK7BEXy0tmtRmkaXMuWnvtv7VpvapUyenTImh0qHt7ivdnZ+ta2y4vs3atd2WLvKiodO3p9//JLiIXn86koWYsoMVHFw1tqZMOGg2o7TZ6snsi0NMlbkNHWXMrIUOPJFtHIm19lfYaLF2tXmzfbc0eO1M7i4iQnR3MwOne2FmJ8vFqWffuK+P2yYYPqo1fjcOpU7Xv5cvH59NG0b2/zMFJTtaNPPw1JoXQ4HEXN8RCscOB9YBd5y40EgJ3oYo/hhb3GPxzHMCDVfj+Pf5500QdI5P8nXbxMESddpKSoOHiW1ObNallYlZo16yDDStPdbVn0QEB1ZtEieywrS9/oNrsiNlZdi54gbtyoL/j0dBHRUNL773sapm98G5vy+TTBwxZ3V3XyKtiqdrVsGeIiXLlSd1jBnDhRKzzFxdmB9u+vvj9rbX35pd5XTo6ooLzzjmYF+v2eATV2rO17yRJ9SAsWSCCg+1u3ttU8fD6NqbVvL7J7t6SlqQX50Uf22jExqvpff+3mbzkcxYDjndYeYcUi4mj2+w+uW9W6JqeE7FsCLOPgtPZGVoSeCNnX1O67P1+fs4GNhRhTwQUrJcUTpXXrQgraimjZi2+/9X6OGhXihhPRbDybjRAIqA54SQ/BWoBWSYITgb2qF9HRmgVo41ZB48kzQhYs0MCQfbn/+aeKWna2vVi/fqo2ubni86lh5c2tyszUG7FrWyUmqo54KylHRakVaH2Z69apflqPon4JmX08bpx6DTdvFr3AoEEqmrt2SXa2WojvvGNvJTVVx/XJJyIJCXLggD6XLl2slbljh/4h0K1biK/U4XAcb07GicODgU+sy7EB8D9gB5BBSJIFOgHYBwy3QvUMsA+tfBEe0s4A09EY2ItAQ2CAtRLvLcQ4Cy5Ys2drQMcyYcJB3j71qQ0f7v386acQi0NElSJEtD7/PKS2oN+vb3K7YnF2tmYI2pwOFaM33vCSExISVLS8Mk/R0Wqabd4sInllmexKJtowxJpaulR1xmbC60Bat/bqP40bpxoYE2MHO2KECk9srDfd6p13rEWUm6vBvfffF0lIkKws1e6PP7bZ/Wlp6ur7/HOR5GRJTNRjX3xhhSmfRbV7t/b93XchRXU7dFC1c8LlcBx3jsUCjt1tzCoZrSSxD00r/y9Q4Uj7K8D1OwAr0cSGXCDGugOvOETbO4FFaAZhLDrn6oxDtKtg7yvGtl3O30xM/gfjLJxLsFevkElXGgoK+anuOU+F9OekSSHHv/32oIoPX38dUlA2ENC3uK2EHgioAeJZO8F6gHaHz6fhHm8RxpwctVasaGZn689hw+zxzEw1YYYM8arBf/hhyMLBqakqOgMGiAQCkpKiP/v3t9ZYQoImZPTpI+L3S3y8xr569w5JgX/3Xd2RmysxMSo8PXpYN+KePbrjhx+8NPi331ZrLytLVEw7dvQsqq1b9ef331uP6ObNuuPbb0NmXDscjmPNsRCsgH2pbwc2WiEInSz8xJH2eTJuhRasYGn17du9n+++G7J4YyCgKuIFqVSUpk4N6eO77w4yvXr29DLRlQED1DyzDB7shamUkHqAIprX8cEHIUuejB2rwmCtkQkT1DjyFgaeOvWgzJFZs0LmVomor/L11724Wj6vn8anWrb0qmSsWKHNZ860x1eu1A7HjRMJBGTtWr2cl8q+fr2abz//LOLzyYYNak327m3vYfNmtai6dxfJyJBt21TYuna1Ftn27Xp/X3wRkhrpcDiOFcdCsGodIkHhX0BrIMq64Qpcg+9k2Y5G0oVkZ+tcKJvOF1wh2Ks4FAio6REyt6hz54MML307h6jUyJH5RGniRPWbWVGaN0/DXDb34qB6gCLq+mve3PP46Yzkli29NbJiYw8quq7q1aGDZ56lpWnSg5fCnpWlmRBduohkZUlurhpGH31k55r5/Rob69BBZN8+z2vYpk1INuH48epmtH7NGTPywmGBgKjgvfGGloDy+2XNGhWuvn2tcG3adJBFtWuXWnydOtlKHdHR+sfBhx+G3LjD4TjaHO+ki9Otu27R0ez3RNwKI1gHVRI6cOCgRRgzM/Nl9/n9atbYRa6Cc3q99G4RtZR+/NH7OWOGvpC97MM1a0KWAtbw0quvhiRjBOsB2rhXZqZaWl56erB67Xffeetu/fKLtvGS74LVKGzsKlg30EuP37JFVcbG3qKjVaM8N2Fysvodv/5aJCNDMjPVEnrvPbvkic+nWZJt24ps2iSBgHo8W7cOSTqZPfsg4Vq9Wh9djx42VX7HDjVjO3cWiYuTmBgtHP/ee9bQTU7We2zfPiTo53A4jhbHPekCaAlkHO1+T7StMII1bVpIiSURzWx46y2vXHtwvUav2oXPp2/ekBV4e/QIERQRFQKv8J/q22uvhbjv4uJCJkHlLV3lCV9wdeBvv/X6GDdORcWbzxVc2t66+PbuVX3w3JSpqWql/PSTiN8vubnqjfv4Y9tHIKBxszZtPFfookU6TluiUJ9Fu3bqzvT5JD4+L7EiOdkOvHt3z53q9+tzeOONkInYc+bojqFDRXw+Wb9e7+Obb2wfQYvqgw9Edu6UlBQNd7Vta5czyclR/2mbNip+Xu6/w+EoDMdUsNC6gU8D11u3YF1gFrCjMP2eDFthXYJ9++ZbJWP5crUwLHFx+UTL79cAzLJlXptBg/IJ37JlGtOxtZmCfXgZfrm56ovzsi/0nf7JJyGL+i5frgpiJwrHxOi734srZWerevzwg2dtDR2q+mG9inmpg3bRxehoNVoGD7YuvIwMdRF++qlIcvJBbkBvva0lS9R8GjNGJBDwMv66drXuzLQ0Fdd33/VWOR4+XMc6fbq9zsKF+odAnz4iWVmya5dq1Kef2tsLWlTt2oksXiy5uapPbdpogklOdkAVrF07dWt6Zq/D4SgIx1qwypM3WdgPpAObgHsK0+/JsB2NGNZnn+VbSHfatIPWGAkKjmclBQL6xp0/32szdqzGYryYVVTUQUsHZ2frO93LIBTReNPHH3sqtWmTnuK5CFNTVR3ssiDBakkffBAS+woKm/X5xcernnqZgj6fxtfee89LaJg3Tz2T1vOoqtGxo5cNmJWlls4779jEjUBAA1WtWnkplNu2qbXUrZsdS2qqClfHjiJbtojfr7UWW7dWl6HfLzrGDh3U3ZiYKPHxqpcdO1qjNTdXFerNN1X1cnJk8WI9pXNnm5K/b58KdceOmj3iCu06HEfM8ah0cSlaFWK8zRJ8u7B9ngzb0RAsv1/fkZ5QiKiy2IK1Ivqub95cSyaJiL4ou3RRy8OydKkaEl5sLDVV39gh1tjPP+tp3iKRmzdrsMwmGWRmqvHlLfshoiZg27bezOJgFXXPBZibq8rx+edeJmFQlLw8kdhYVcyePb2ahME5wF5+w6pVahrZ9L/kZNWGjz8OWb5k9GhtY7M9glOqvMz09HR1FbZrJ7JypQQCqv9t2micLDvb3sBHH6nybtsmWVnqeXzzTRU5n09UTTt00IcVEyP79+vX4MoqgewcrWPYtq3ek/fXhMPh+Dv+TrCM6Mv1H2OMqSAiKYc5VhsYA/wkIh8fUccnGcaYpIiIiIikpKQjP9nngxIlAMjJgbfegvbt4dxz7fFx42DXLmjeHIDERHj3Xd3OPtu2GTAAAgF48UVAm3/5Jbz/PlSpgh7r2hWqVoWnngJg1Sro10/7qVwZyMqCzp2hRg147DEApk+H8eOhQweoVAlISIAuXeCGG+DBBxGBX3+F5cuhbVuoWBHYuRO6dYMGDeDee/H59Dp79kDr1nDmmcDq1brzjjvgzjvJyDT06AHp6dCypW0zaxaMGgVNm0KTJsTGGXr1grAw+N//4MwzBP74A2bP1ja33cb2HYb+/eG00+CVV+DM03JgyBBYswbuuQcaNGDVasOwYXo/L7wAZ5ZMhUGD9KHdcw9y8y3MnmMYNw7OOQeefRbOzImBn3+G+Hi45x78desxfoJhzhw4/3x9pGcmboNhw/QmmjaFm24CY47834PDcYoQGRlJcnJysohEHup4QQQrGugIDBKRwCGOtwQ6iMh5BRnwyUKhBGvqVNi7F557DlDdaNtWhaRKFdvmzz/1hdqyJQAZGSoir78O1arZNuPGqQp16ABhYaSlwXvvwdNPw7XX2jYTJsCCBfDOO1C6NCkp8MkncNddcNttIW1mzYKOHSEigqQk1bHrr4eHHgppM3WqDrRqVRIS4Jtv4JJL9DbCwmybyZOhRQuoVo34ePjuOzjrLPjPf6BUSclr89xzULs2sbHw449Qvjz8979w+mmiijllCjzwADRowIEDeML1n/9A5Uqi9z5tGjRsCPfcQ8yBMPr00b8FXngBLjw/oP3MnAm1a8Njj7EvtiQDBkBmJjz5JNS8zK9tZs+G6tWhWTP2JJbn558hNRXuvx9uuMaHGT8O5s3zM7q8IQAAIABJREFUlGp32pkMGQLJydCkCTSol0vYxPH6nCtV0s69vz4cDkeQYyFYC9BK5tuBnujSHJvssRJo1Ym7RKR8oUZ+glMowQIYOFDNgocfBvSP9Hbt4OOPrdUCKhBLlqggGUNururO44+HCNKqVWptffIJnHaaZ1idfTY0a2bb7NgBX30Fb7wB1aohAr/8AlFRat2VKgXExqpK3XGHvoWBiRPzNKpKFfQt/s03+uOVVyA8nMWL1RB54QW45hogOxt69ICUFFXXyEg2boTevdUAefhhMH6fnrRhA7z8MlSvTlSUitJZZ+mucmUFfv9dhfSee6BRI2LjDD/9pFbpyy/DeeeKCtK4cXDFFfDkk6Rkl2bAANi3Dx55BK67Dli2DEaO1Af73HNknlaZoUNh3Tq45RbtvsS2TTB0KIjAY4+RW/3f/PEHLF6sovzEExCRvFutt6QkaNwYf/3bmDo9jBkzIDJS/7tcXP6A9hMdDbVqqeiWP6X/V3E4PI6FYBl00cQP0KKzAmSjZY0qA+XQIrR3FmbgJzqFFiyAH36ASy+F228H9B3fsaNaSVWr2jaLFsHYsfDRRxAWRiCgHrrq1eHBB22bAwdU6d56Cy66CFBDZu5ctdrKlgVyc+GLL9Q8e/JJALZuVQvoP//RdysAI0ao+65tW4iIICVFNeq889T7GBYGrFwJffvCM8/ADTfg96tmbtumGlW1qh1Tt27qY3vpJShVipkzYfRouO8+aNQINS379lV1eflluPhitm3TXVWq6K7y5UKsqdtug6ZNSUoJo18/1dhnnoGaNdExDx2qFs7zz5NboSKjRqne16mjz6pEXIwONCVF3ZvXXc/8+fp4zzhD+zo7IkOfwfr1KoKPPMK26HKe5+/OO+GWm/yY6dNgxgyoUAEee4yEMy5hxAj1jl56qYplxK7VKrpZWXDrrdC4secKdjhORY66YHknGlMaeAItQlsbXawxCy0k20JE9hZsyCcHR0WwQNXn+us1/kOe6+/NN+HCC22btWv1Lf7ZZ1CuHKChk5gYaNXKhk1yctTKqldP36qoBdW5s4pI9eq2r0mT1Gx6+2044wz8fujeXQ2Lli3t+/TAAfj6a6hb11PFpUs17PPSS3DVVWiM7OefYeNGvcDZZ5OcDN9/rwbF//5nh7p+vcaurrsOHnsMMWGMGaOG0WOPwY03okrQp49e1wrXjh16y5GRuisyQlQgxo2Dq6+GJ54gy1+SX37RIQRdnCZ6H/Tvr36/Zs2gZk0WLdLQ1xlnaHyqyhk5KiRLlqjaPfYYsRnl+eUXfab166uhGb5+Dfz2m97rfffhu+paJk02zJ2r43riCbgwIkmDelu3qin26KNsjj2D335Tl2HdunDX7X5KL5qtohsWpn+g3HSTVX+H49ThmAmW468pjGCJhMTmRVQcrrzSs7Sys1W0/vc/uOwy227fPhWsEPNr4UJ9n374YYjXaehQ2L1bra3wcHJzVRPPOUfDRsagiQSdOqmwNW4MqHusZ08VpKuvtn1NmqRb69ZwwQX4/SoiByVTpKSomVa+PLz6KpQty44d2leNGmq1hIcD8+er5dKoEdxzDwEx/PqrGpCPP65WkCdcMTE62MsvZ+9evaYxOrZzzkEzPkaMUDPshRcIVIhk4kRNGKlZUw3IspKhz2LDBhXxe+/lQEIJBg2CuDgVpAYNwGxYr335/fDww8hVVzNrlt72aafp2KqdlwVjxqhr8Zxz4PHHiS9ZleHD9VFXqwaPPgoRcdvU9ZiYCNdeizS9h0WryzJxohq49evDbTfnUGLGFI2JlSqlA6lTx4mX45TACVYRURjBmjJF3VlejAnUNLnkEs02Q5MHPvpI32c332zbpKdrEOvll9VdhYZKPvtMLa1LL7Xt1q/XOFL79urLQ/MKfv9dDatKlWy7ESM0Bta2LURG4vdrrCkhQcNd5cqhJt/330PJkmqClS5NfDx8+63mFbz0kh5i+3ZVqSuu0KyPsDCWLdNYWb16mrwRZkStjLFj1SS6/Xb8AcOIEWrB3X+/es7IytIY37ZtqgTXX09CghpO8fH63K64AjUhBwxQNXjqKbjsMtatU+uzVCkVy4suFDy/X2QkPP00/qrnMmmShseCsb6zKmRphuKqVWraPvkkSeYMRozQYVSvbt186fv0uUVHq6l5//1s3lue335T7b7uOmh6t1Bm3TK1BrOzoV49/Lc1YdaCUsyYoX+j1K8PDepmUXLGZP3Lo0QJFfObb7YK73CcfDjBKiIK6xIcPFhfXE8/HbKzZ081W2yKuYiGgc4+W9/bgFoCXbp47ifQ93WnThqH8uJaGRm6s3ZtL9UvOVnDWHXqaC4AoMrZpYuqyv33A2o1fP+9JiTcd5+1yrZu1XQ+m7qOMZ6nsmFDb5e62YYO9bL3MIbZs1ULGjdWPTaImjCTJ3vtAmIYM0aFtWFDbRcW8KnFsnix+vvuvpvs3DCGDlUvaf36tl1muiZDbN6sO++6i6TUcAYP1kTLW29VfQxPjFMFjY7Wm7vzTvbuL8HQoeqNDLoCS0TtgOHDNbni5pvhjjvYtL0kI0dq3ondRcn1q9TyysiAm29GGjdh6epSjBunHtp69aBJowClls7Xe/X5oF49fA0aM2thaWbO1P+cN90ETernUHredA08gvpKGzWCMmUK9O/L4SiOOMEqIo5GDGvoUBWbZ58N2TlihL49W7Tw/Ia//qoi8sYbIZ6j33//P/beNDqu67wS3bcmoFCYCoUZKKAwzwBJgDMlTqLFSZGHSFbHzui8jpNO5/XqOJ1ppfPiTtzr9XN+vLys1U4nK467O5HjIbJjSZREURzEESQIkpjnoVCoKhTmscZ73o99pwLpJLaU0FbwrYVFqe695557btXZZ3/f/r7DXKPf/m0tkP+d79Br9Tu/owgtADKLGzd4XmYmACq5L1+mKjEvTzlPFTYo7j+AKsG33mI6WFIM7MIF4Od+TvMdvvMO23zpJU7SAHQmde4ccPIkBCRcusTbHD1KIDRJhrhURwfw0ksQZguuXOGlLS1076XYlPMuXCBQf/azEI50XL3K+5aWkk05swUR78IFPthnPwuRX4Br1/gcWVlkU2VuAbz/Pj9MSwNeeglyTR3ef5/PYrMR41uaZI7d22+T9Zw7B7l9L27clPDOO3zMU6eAI4dkmG5eZ2wwkQCOHYN89Dhudlpw8SJx6uBB4NSJBFLu3+J50SjQ3o7EqdO49ciBS5dILFtbgbPPJ5A1cJvnRSIc/HPnDC9rx3bsx9N2AOsp2YcluvjGN+iC+/znDR9eucK/3/s9zT2kute++EUgI0M5b2iIFOy3fktz/c3MUGiRpPwLhZhV/NxzpAYg2/rylxl/+ZmfUbBxc5PxKLudKGWzIRKh3HxpidoKpxOcgb/2NbKuX/olwOOBLEOLSf30T5PYQRgUfmfOcHaXJFy9SmLS3k4yabGAzOxb36LK8Wd+BnA40NND915WFvExPx+85//+33yuz3wGqKmB10v9x/o6GWZHByAFA6Sxc3NkZ6dOYWHJpIX49uwhy0yNrLDjw8NaMGrD5sRrr5HF5eezj6W5YaLj3bt8AR//OKLVjXj3XZIiq5VDe6AjDtP7V/n+hACOHIF87ARuddlw8SLxZ/du4MxpgYzhLoLm+jpQUwNx7jx65gpw4QLfT2EhSWplbIj3DoU4GM8/T3fkTpLyjv2Y2Q5gPSX70FSC4Kr+3j2SIG0O6usD/vzPiVAKMwqFGNdKUv1tbQFf+hIrUbzwAgBoyj9ZBv79vzcoqb/9bd7oN35DUUww1/XVVwmYjY3KecPDRKmjRzVfXyhEN6GagpWSAsbUvvIVBm9+5VeAggLE48ST3l66O3ftAifuCxcYvNPolQn37hGwKyo0jNIByWIhSpWWYm6OoaqlJQLSvn3gJP83f8PzDx0Czp9HVLZouVOVlWRTWRky2dm773IcP/1piIpKdHdTOZhIAGfPkgFJY6O6aOLAAeDMGQSWU/HNb1Lz4vEwjyzXtsqL+/qI4C++iLCnHu8o4Sirld68wwcSMN98n/ePxSjEeP40How4cOECH8Hj4RAXrY8Q3INBvpvTpzGb04w3L0gYHyfrO3oUONy4BNvVi4y1AUS/U6cIZDu2Yz/itgNYT8k+cKWLhQVK0BTr6mLY5I/+SBExAJSz/cEfEHUUhIrFSJYqK7V0KtprrzE/6jd/U5O+q8q/X/kVKvYA8L5//MdEp898BpAkxOM8b2GB4o1s9auk+voM7r+hIQr5WlsJCGYzSBG/8hUyr1/+ZSAvD7GYXiHppZcUFaBQXHZ///dkCK+8AthsGBkhQ7LZmIBcUgKCxte+RsXguXPAkSOIxSW89lqSGp15Wrdu0YeYns5BUfK5vv51gsKpUwy9mVaWOMgTE5qKImrPwltvsYnMTLoC62oFUe/CBYO87wQmvBZ8+9scp8pKgmeueYngNTTEBs6dQ6S2Be9dlnDjBofx0CHgxHFFiPH22wT6igrg/HlMhIvw+ut8zPR0CkV3ly/CdPFtor7JBOzbh8jhE7jW5cCNG+xSeTnw/CkZ5UsP+H1aWWEDJ06Quu7ke+3Yj6DtANZTsg/MsP72bzkpG3yBY2NU3/3+7xuUfPE48Id/yFwtRUEIQFO5/e7vGiTts7NEs5df5iypXP6nf8pJ7td+TWFGAGM43/oWXXoKtZqbo0ewooLAYTYrDfzP/0nW9fnPa4nJnZ0EhAMHqJ4zmUCA/bM/I7X7xV8EiooQN+gmzpyhV1KSQHD92781FPjL0ZSAc3PUfxw8CEhygszj+nXGr37qp4CMDPQranRZJutpawPH89VXmb3b1gZ88pNI2Oy4eJE4mZ7OvtbWggDz7W9TRXHgAHD6NJa3UvDaa3zU/HyCV3lpggN9+TJvdvw4cPQoRqes+M53CF5lZXQvFqWtEOAfPaJY4uRJJPYdxM07Zrz3Ht2B9fVkdLmr43wuv58uxo99DKtVu/HOuyY8eEBs37WLca/skbvU7G9ssFrHqVOYTG/G2+9ImJ4mNu3dCxzds4aMe5fJohMJIv/Jk3zgHffhjv0I2A5gPSX7UFyC77zDnKLf/E1tQlldZQHbz33OEIMCODsPDTFepVAwv5+qv5/6KYXBAJzpvvY1Sr5//dc1tjU+TpeeUp6PFovR7RgMUmzhdAIglnzta1ysK0I/xre+8hWu5D//ea0K740bBKRDhzjBa4zrL/6CYPALvwBUVEAIvdTTrl0klzYbWFPRWOCvqQnxOEnYrVuc4F95xeAufPVVzvwf/zjQ3o6tsKSp0Ssr2a7TqTzEd77DZzx9GjhyBCurEr71LQJSaSmZX2GBYODtwgWC87FjwPHjCC5Y8NprVBkWFJBNlZfEGZu6epXgdegQcPIkpoKp+O53yZJcLo5ZXXmYsbtbt/hO2tshnjuFodkMvPkmsT07m+GoVs8qK2d0dfHctjbIJ57DQ28OLl5MJk8dnnmY37tI+gwA9fWIH3sOd72FuHKFjNLhoAhyf8kMbNff4/cGYIzu+HFt0bFjO/YvbTuA9ZTsQ4thqYGcP/gDTdonyyRK1dVkBJqNjrJQ4G/8hjbpCEFsmJ8nPtlsyrleL899/nlNaCEEb9XVRbZVqpYvDoVIrYqKGKBSGnmi+m9lhSwqHCY7U6r13rhB6XpbG3HHagVnz69+lUqQl15SCvvp5Conh/nB+fkgYH3968wh27+fFMtqxeAgP05iUpEI3XBdXaQ3r7wCuFwYH9fV6M88o0jPESNS3rzJWf8TnwAaGzEzQ6ANBkncPvEJwJUV19kUoCRLHUNggWxqaop9fuEFoL4mwTbfe49j0dwMnDmDeTkHb7xBULTZiH+HD8qwPLrPAV1f55idPYul3BpcVMJRkkQgP3lCwOl9RGRfXCT6HD2K1fp9uHzdiq4ukqeyMuC5kwKV0UFIl97lg5jNQHs71vc8i+u92bhzh3idlQU8+4zAnqwxWK9f5oMIwS/YsWP8Lu0wsB37F7AdwHpK9mGKLuDzUTjxu7+rlHKgfe97jNf8zu8Y0nEiEaJZWZlB3qe7Ez/7WQPbAhjbunmTDEqpIL6xQXGh1cr4liaBHxgg4zL4+WSZE/vt2wSivXuVcxcXee7WFumg2w0A6O4mwJSU0NOXkQEyl29+k+B8+DBFFxYLgkEyuaUl4tP+/UqOVmcnASk9nfTR48HWlp7XW1VFJpWdDU6+3/gGfXMHDwJnzkC22PD++9R4mEx0we3fD0jra2Rd/f1Enk9+EqiqwsgIh2lxkcnXL74I5GYr4HXlCtFy/34qDTfteOMNloOy2+ni3LdXwDzYR3RfWiIgnTmDcFktrlyVcPMmh6C+niBaIAJkdGNjWoxKPnocD0cduHiRgJueTrzc37QOy81rZIHxOKWDJ09iKq0Bl96jIANg28efiaMk0EX/58oKX3BHB5abj+D9R1no6mITmZnAoYMCe11jSLl1lTE9gC/tyBGgqWmn8saO/bPYDmA9JfsggHXpEiclpVA7bWuLMvZPfUopsEebmmJe73/4D4ZtRQBOSq+9RhehwnJkmd61qSlWZtLk7+vr9AempVEUoQSy1BJKjY0EOq3AwvvvM76jZfpKSCTojevuJpZpXVxdZfbw/DxFHEo8bGKCYGQ2E1fLy8FV/a1bBKPCQurfc3MRi/GjO3cYbnnlFaXvi4tEv8lJsrOf+AkgNRWjo8S/lRVO6qdOARazUtHirbfYr+efBw4dQiRmwptvEgPT0siO2toAaXGBCDg2Rj/exz8O1NRgZIR9WVggA33xRaC0WGbn3n2X76muDjh3Dptpubh0iW0LoQv2MreC7MfQEANMBw9CHD2GoWk73n6bMbrUVPb9QEcctgedBMaNDdKhEyewXrMbV6+zGn48zld84gTQlOOHdPk9Li6EoLv12HEMRitx+YqE2Vk+fkMDcPRgFKXBLr7PlRW+jF27sNL2LG4O5+LuXaaDpaRwIXKobAaZPTfobhSCqLl/Pw/uVJzfsQ/BdgDrKdkHZVjf/Cbn4S98YVtdQVVt96u/qh2IxVgxPSeHWgbt/LU1sq2aGk7+yoFAgN7A3bvJRLTzR0bYvhpwUg7cvcuUJaXMHz8WghP0G29omy6qwPXtb5NxJYkowmE2MjBAkDt2DJAkLC9TszEzw7afeUY53+ejfH15mQeUzQ+Hh+nWC4f58YEDCuu6f5+BLVnmgX37IAvmdF26RFw4d07JwYpFqUq5eZPIcO4cXWUbEl5/nSzN4eDHu3YB0sI82x4ZIVKePQu0tcE7I+F732Pfs7L4vC0tgDQ8xHGZn2fA7PnnITe1oPuBhIsXieG5uVT8NdXGIN25Tba2uUn/56lT2PQ04v3rEm7dImgUFXH8G4qWIV25zOeVZbKe48cRdNbjvcsS+vv5cVkZAa/ONsHzVapVXg7xzLMYSNTi6jUJPh9fZWUl8MzBOGo2H0K6/j77DgDV1Qh3HEHXchVu3pKwqmzdWlMDHGpZQ9VCJ6S7nQRTSSJYHzhAmrvjRtyxH9B2AOsp2YfhErx3j6zli1/ctoDt7KQm/D//Zy1fCiBIvPoqNRoGzyEDSN/4BuXvBhp27RrdeVqFddUuX6Zr7NOf1tSEwlDm7+xZLceXBy5eJGN49lkth8ooomhvZ4jKagVn0wsXyBgaGpRKtHYkDGK/igqSscxMEJzfUDZILCmhC1CRxb/+OolNbi5ZV2kp6BJ94w2OkdNJuldVhXCYH3d1cSxffFGpN7i5yQP379OHd+YM0NGB9Q0JFy6QMdpsBItDhwDzxiof7OFDouAzzwDPPouVLRveeosCQJOJc/bx40BaeJGxqZ4eDlh7O3DyJOajmbho0Ea0tPAeufIcB80gmsDJk/BLxXjPQJyqqth+hc3H9zU4yANlZcCxY/DaCUhDQ/y4qIhdbU6fhOnG+9AO5OVBHD6CiaxduH7bgpERfux0AgcPCLQ7x5By7ya0A5mZEPv2YyR7L249TMPoqE62OvbI2O8cRkbfbTJTgAuCPXvIwjRp647t2JNtB7Cekn1YMSy1eO2/+3eGXCmAsZAvfYnIoVRxB7jQ/fKXObdr+1MBnMiN+nVFHRiPU/cwMcF7aBvhyjKpzL17ZGdKnpUQDMW8+y5vfeaMYSF99SoDa01NBBbFtXjvHllXXh5FFNoGlP39RFirlecrYDo+Tjze2jKwKAmkMn/zN/THHT7Mm1utmJujZ3B2lgv8T35SyZNdWuKNx8boM/vJnwRKS7G2puf1ZmTwHi0tgLS1STC6d49gdOIEcOQIIgkLLl0iIROCLO3UKSA9JUZ32rVrpEE1NcDZs0jkFuD2beLI1hZvffo0UFMlQ+q+T+RfXeXNjx+HaO9A74BZS7+z2ehSPXxIIM07xMH2+/ky29ogjh7D+GouLl8mC1f1EUePAhXmaUjXrurKv4IC4Nln4Xc14/2bZvT0aLiDgweBjrI52LtvUukSjxNg2tuxWHsAt4ec6OriV8dsJsAfaFhB6WwnpK57BHsAqKzEevMBdK3X4W6XCSsr/LiwENjfsolW+QFsD++StUkSb97eTiDbSWjeMYPtANZTsg9TdBGPM05VXr6tgjtAUOntZazKQMPu3qWrLSkpGKA68E/+hCvel17S0GZtjftFyjKv0ZKD48ruv/39BK7WVgA643rzTQLKJz9pyEXt7aX7T9nBV61xNzvLuNX6Os/XdkVeW9OrUuzdS/qTkoJYjOTn1i3Ou6+8ojBHocSjLlzgJH7uHMtbSBIGBxl6Wl1l++fOKdjs9xO8fD7OpJ/4BFBWhtVVMrWeHs7Vp07xeUzxKBHn+nXK7nbvpmsvPRNdXSRNGxtkLmfPkvFgZIQDEgySrR07Bhw4AP+8Fe+8Q2WgJJHNnjgBuKyrvIcq7SsuBk6cQKSiHrduM7F4a4v9P3wY2N+RQOrQQy4M5ueJIq2tEEeewdh6Aa5u00c88wzQmBMgo+rt5T2UuNNK3T7ceuRAVxfvoQLSwd1hxrXu3CHgCwG43Uh07EcvmnGnywKvl/fIyADa9wjsc40hY6CTIJlIcAHS2Ah/2X7cDbjR0yshHOY1bjewr3YZTZH7sPZ2Q0O3jAyO8e7dhhXNjv1rsx3Aekr2gQDrvfc4iSvV0VV7+20e+u3fNgAKwEn4j/+Yy3gD24rFuIvI8jJFGUoFJ9q1a5zAX3rJsD8J5/WvfIXn/tIvcX4DQAbx139NWvLyy0r9I9rt2wQJj4f4pOGm30+EWl3lfXbv1pr6u7/jPN3YyEPp6eDk2NVF+mM2E1QUX2UgQBbl9xMztZ3lIxFoqon0dM3XJwTdeW+8wQl5zx4SModD6ddrrxG8nU5SrIYGbIUlvPsu52ohiJ2nTinVMh484AtYXSUAnzkD1NVh1i9pgj6zmWB37BjgkDbp9lTVe0qSrlxdi4ePJLz3HjUjNhuvOXIEcCxvc+95PMDRo9gorsHNWxJu3+bjpqbymoP7EnCM95DlBYMc8+pq4MgRzKRU4cZNCX19XIRo+oj6NaQP3OV4qXGn2lok9h5A71YVbt+RMDPDprKyOAYdhTNw9N4hqicSXCQ0NGC1fh/uL1Xg7j09tuVyAe1tcbTb+5HWf09HUIsForEJM0V70RUqQ2+fhGiUhwoKgI7aVbQkHiBt8D4HBuDgNDXxO1BRsRMT+1dgHznAkiTpJICfBnAQgBvAIoBOAL8vhOgxnHcFwNEnNPG3QohXtrWZDuBLAF4CkA2gD8AXhRB//wH6+cEY1t/9nc6ctOQpLqz/63/lqv7kScP5QhCAOjuZh2Wo3B0IMI2qvn6b2k8IvSLtz/1cUiby9DRV6U4nY1ya5yYe15O1nntOE1sAJEj/63/pyj8t/zQSYbDswQMu419+WdPKDwywC9EoFXoKUaK76bXXGCsqLeU1ysaUjx4R0zY3GVf62McU7+PqKsURfX3s8PnzQFMTBCR0dxPXNjbIOM+fV8J/i4ukWAMDbOTYMRakNVlw9y7Dc5ubJD9nz1KcgLk5XeVnMpHKnTiBeFombhv0E9nZbG7PHsDsn+FqY3iY415XBxw/jki+G3fuMES3sUEw2r+fz+UITbKx4WGOY3Ex8Oyz2Kxowp27Jty+TbZqsRDEDx8SKNwYY2Ojo7wmKws4eBBrNXtwt9eOO3d0T15NDbB/r4wajMB055Yed7Lbgd27sVy7D3fHXbh/n/cBiLsduxNotQ7A9uge/bdCaKwq5NmL+/NlePBQ0q5xOoE9LTG0p/YhY7SbIKbuUlpdjUBpBx5s1uJhnwXr6/w4LQ1oq49gT2o/8mcfQJqa5DVC8PvQ2srvkiZ13bGPgn0UAeubAFwAvgFgAEABgP8EoBnAMSHEbeW8KwCKAfzMtibmhRCj29q8CGCP0s4EgJ8D8BkALwgh3vwh+/nBXYLj45T//cf/qMyUNCH0HNovfGFbGGB5mRLAggKWZDfUjFPVfmfOaLnCtFhML6/08z9PZFPM66UqPSWFwJWfb+jEpUucuOvr6atU4mKrq2xueprCgI99zACSPT1EKCHoFzSwru99j3ibl0e9h5K6xdjVN79J5G1oIPPKyoIs01349tvaVlJ47jkFvFZWCES9vezXxz5GumAyYXCQhxYX+TznzimbW0YiZEU3buiJUadPA/n58Pn4qOPjfJb9+xkzSrcnCN6XL5MVOxwUn+zbh+UNK65e1T1+hYUcj6ZGAWlkmNeodKauDjh6FFt5Zbhzhx7PzU2+vvZ2PltOeJZsSpWVZ2QABw4gvqsDj0bTcOMGsRSgq/LAAaDVvQTLvdsUlYTDGqMS+/ZjRFTjTqek6SlsNpKZvU2bKAo+4MtQ405ZWUB7O2aL2nFvOBOPHvGdCaGAWFsMLeZ+pPR26YE1kwmorcViRTsebNSg+5EZ6s/BbgeaGmS0Z4+hJPQA0tAgNMrldGKjZhd6zLvwaCYHMzM6vpWWCLQyLZPeAAAgAElEQVQX+tAQ70HaWI+OpGYzX2JTE8dyZ5+wH0v7KAJWvhBibttn2SDQvCeE+JTy2RUA2UKIXY+3knTtWQBvAPikEOI15TMJwPsAXEKIhn/o+n+g3Q8nhhWNErRyc4kYBrdIKESBxb59SSp0Wk8PKdKLLyZRMaNo4uWXk1K6OGF/9auclT/zmSTp4NwcD62tkT1p1eABXTyhbePr0e51+TIBJS+PgkBN1BEOk0F1d3PF/MorGhrOzTE0NzNDD9enPmUQQw4MkF6trHCF/cILQGYmEkphiUuXtK2k8Pzziqtxc5NUqbOTg3TgAMfEbkcwSLfh6Ci7f+QIYz8pNkEG9dZbHGgtMeoA4iYbOjtJfjY2iBtHj1KMYQmvE1TUMhIuFw/u2oVAyIwrV3S8KSjgoeYmAdPIEBv0enmwvBw4cgTRynrcf2DCzZsEWIWU4OBBoKZgFVLnHYpENjd5sL4eOHAAsykVuH1HQk8PsVeNUe3rkOHeGqYUXZX42WxAaysiLR14NF+Mu/ck+P0cbrudX4OO6mUUzHTxfa2s8F4ZGRCtbfAVtuP+lAs9PTomOp3ArpYEdjuGkT31kGMZj/Ngfj626nahz9KGh+MZWmENSeL3Y1f5EppFDxyjDzn2ACBJEGXlmMnfg55oHfrG7Vhb03Gx2hPHnowRVEf6kDIxyO8yQOZXU8OFTl2dIQt+x34U7SMHWN/PJEnqBCCEEPuV/7+Cfxpg/TnoCswRQsiGz/8PAP8DQJMQov+H6M+HV+kC0DXrX/iCgXrQ3nmH+oPPf56/Sc2E4OR+9SpdfgYAUven6uwkVmgVKgAklVJ/4QUyBgUNNzeZHjUywjn/Yx8zKBGXlvDYNr4Kw5ubY/dnZ8lOzp0zFNr1eolQoZC2rbwaCBsdpadzeZlzzosvKoxSCLKn118npaupYV/z8iCUtCx1Kyml8DkFG4kEx/LSJU5qJSXsZ0UFolGSq/ff56GiIj5fTY2iILx2jdfGYkSb48eB5masrnHX5Hv3OCerGNXWBpiXF3jdgwcc9OxsouKePfDPW3HtGrQ4k+K9Q0cHkBqcouBDjWdlZAD79kF07MVYMB23bunexfR0vr+97TIy/UPso5p3lZYG7NmD+K4O9AVcuHMHmmgiNVUBo5YICuce8QHUzOKUFKClBZsN7Xg4X4L73RICAX4NrFYSmT3Vqyhffkjlo8rErFagoQEL7l14uFqBh71mLC2xSYuF38/dxUFUrT+Eub+H4CcEY1xV1fAV7EHPVjV6h6xYW+N1JhNQXSljt2sa1eFepIz26X5NSUKivBLjmbvQF63B4JQdm5s6kFWUxrA7fQTVsQE4ZoYY0JQk/pWVEeDr6nYk9z8i9q8CsCRJygMwBeBVIcTnlM+uANgPIAwgA2RgXwPwfwshYoZrb4FAd2hbm/sB3AbwaSHEN36IPv3QgKUWKT1/ftuBzU3WV8rMJDoZ3H3RKMUSy8tUrSeJMtQNFUdGqHU30KNEgljR1cX7Kfm8NFkm/bh6lYEYLZlK9wi+8w4J0mc/a2BBQiTv2PvKK5pkXQi6Jt94g82fOaNUXVdzuowBqv376ZJTVsV9fQxRra6yuRdfNMwzw8MEr/n5bX4+zt1vvsm52OGgQm/fPsVN6fUS7ScnOZ4HDhBtHA7MzpKYjYzwFk1NxKjCQlDkcPmyTpdKSzl4dXWYX5Aew6hnnuEQWtcWCUTd3XwvdjvRZv9+LMuZuHULSVLylhaOT3HGGlcXRrFEdTWwfz/Wi2vR1W1CZyc08UNREYGvtWoDqQPdBKPFRfY1KwvYvRtbDXvwcNqJe/d0zUZKCsNDuxvCKF3uhXS/S3dbWixAbS0iDQSHh32WpHBUWRnQ1hBFs2UQ9uGHXG3IMg/m5iLW0IphexseeZ0YHuZ3DyAWN9XF0eYYRdF8D9mm6m9MSUG8up6AtFGOwVGrhlVmM1DlSaAtaxLVsX6kTfRzbNSvb1EJppy7MCDXYdCfhZVVSetrTraMthwv6jGI/KUhmBaVRGlJ4pekpoa/k6qqHffiv6B95AFLcd/9HYAzAHYLIQaUz/8LAB+AQQDpAD4O4BcAfFcI8QnD9cMAhoUQ57e1WwNgGMCvCCH++xPu+48hUVZWVhZ+WIb13e8yPvNbv7UNfAAyiz/7M8aNknx6nHj+9E/pgvu3/3bbby0cBv7yLxlc+tzntAkd4Lzy+usEy4MHGSZK2jKpq4tUJyuLPkGlGjvAOf+v/5qL5dOnOTlrrGtlhYg4NsYf/0svaVXfo1Fi2u3bnLe1skiS0qHOTp4QiXD2PX1aY14jIwSvhQU+6wsvGHKi5+aIUGqFWc3Pl4L1deJMZydvUV1NJWBpKcic7tzRVRMZGQShjg4IswV9fbw2GNRSovDss4q+xevldWpCbmEhDzY3Y3HZpGFULMZ3omAUsm2bBJPbt4k2Bp9fvKIGPb2sdqG66LKzCbZ7dgs4AmPs7/AwH8Ziod+vowN+axnudUl49IivXcXU9nagpWyFINbVBS2o5HAAra0IN+xGT6gQ3Q90taDJxFfX1hRHvXkYtn7FxScrDgmXC6KlFdPONjz05qCvj0QGIADW1wNtxSFUrPfAPNCrV9GQJObFVbSiH43onaJ7ULXMTKCpOoIW2xCKl/pgGh8lyCuMLFFRjYmMVgxEqzA4ZdeGTwhW2W/NnUWdGETe/ACk+ZAu2khLw0JhE4YsTRjadGM6YNPA02QCKvPX0ZwygsrECLJCo5AiYb1TTid/N9XVpO6am2DHPgz71wBYXwbw6wB+XgjxV//IuX8E4HcAPCOEuK58NgxgSAjxwrZzVcD6ZSHEV57Q1j8rYAH8Xf+3/8aA+zaFOyeLV18lI/nVX33MTTgywhBWfT3DStqmjwBnk69+lazilVe49DfY7dsML5WU8FoFX2ihEBUVwSB9gqdOaegUj5NxXbvGaz796W07VYyN0Q+5vMx7nj+vCTU2Nii6ePiQH509y1M05tWlbBe/scEJ49w5rUZiMEiwHRsjPj37LDHKZgMe8/MVF7PPNTWAJGF0lH32+bhi7+ggwcrKAgHk6tVkX98zzwC7dyMhWfDoEZ81FGI/Gxp42O0GEcbo78vI0OrubZkcuHePWKN+PSorSe7q6wQn5lu39BhTaioTt/fuxZI1H3fv4nHlXgfQ2hCju+zePWgzv+K/E7t2w2erQNd9Cb29eognJ4dNt1WtwznTQ1RV0VERMiSaWjFmqcODfhsGBwm6AAGluZlgVLTQC6nnke4alCSgvBzh6mYMWZrQO+HQCBfABUp9nUCr04vy9T6Yh/r1nCwAKC7GirsZg6ZGDAScmJgkOwKIETWeGFoc46iMDCBlalgfDADClYtgfguGpHoMLRfA5zdp11qtQFXRJppTRlARH0ZGYFRHVwCJTCe8Wc0YluowulGEwAJ/OJKkALdzEY22UZTHRuFcmYApGtH7nJ7OF1lVxX93kqJ/YPtIA5YBgP5PIcSf/BPOrwUwBOA/CSH+H+WzHzmX4HZ76y1Oqr/2a0/Yqmhjg5QqkdiW8Ut7+JDsp7aWxMigkOfM8/Wv03f13HNUKRiqcE9PU6a+uckiEYqgjybLlGm/++4T0Wl+Xq+HWFnJ67V8UDVJ6vXXOVl0dNA3qIDX+jq9dPfvc4I5cYKgrYHuiLJdfCjESeHUKaXon4RIhNh0/Tqxyu1mHKqqSrlW9fOpkvSmJt6goADxOHHx6lVildVKJnT4sDL3LCyw8e5ujndGBvXne/dCtqViYICHVWZSVMTDLS2AZesJLr3aWmD/fojqGkxOkUmpxMVq5XV79wIlOVuQHj6gL1UVIaSnE9H37MFsxIWuLr5rVWiXn8/31dYYQ/pUH/ts9N9VVgJtbVgoaMTDwRQ8eEDwVBXq9fVAa1MC1WIElv5H7JiKVA4H0NSElbIW9K260dNn0p5ZvXdTg4yWrGnkBnohDfTrgKL4Djc8TRgyN6J/JhOjo7p70GYDqqsEmnP9qI4OwD41yBWJ2m+HA5HKBoymNGFoqwxD41Yj3iDXJdBUMI960xCKVodhmpnWUVKSECt0YzKjBSNyFYYXcrC4ZNK6ZbUCNXnLaLCMwBMbQcbCJCQDICXs6fBmNmHMVIOxSAl8i2kwTp9FGetoTB1HpRhDwcY4bFsGADaZuFiqqOBfWdm2H+OOAR9hwJIk6YsAfg8G8PknXFMPSuG/IIT4Y+WzvwDwKVARaBRd/CKAP8ePiOgiHOaWH0KwhNJjxbEDAWYJZ2Ux43fbCT09FEskbeuhmhAEn7fe4sz+mc8knRCJ0Buobt3xyivbkpAXF+n2m5zUayMZgHNsjNcvLnIi/PjHDYdV9nThApGxsZHsSQmIRZWCE+rW7y0txFUN/FZWCEDqNrwtLQQw5YTpaR4eH+ek1NrKOFRennLvvj4++9ycXobi6FEgLw/RKMnK9esEMNUNeOSIQu5WV8mE7t7Vg05tbfSpFhZidpaHVaWeIsbD/v1AUaEg8N6+rTMpNYDU0YFobjF6e9m0CgapqTy8Zw9jWtKDbqK6Gpuy2/n8u3cjYC7Bg4cSHj5EUm5TSwvQ2iyjND5JEOzv11HO4QCamxGta8HQJitUjIzoYJKZydfT7FmHe7UPUn8fkiR+ygmhgmb0LZegb8CEYFB30zmdQH2tjObMaZSs9MM0PAhNkSFJgMuFaGU9xlKbMLhahKERsxavUsoeosG9jgbzMIrXhmCaGNP7DkDk5mGhoBHDpjoMrxZh2mfW+g4A+bkymrNnUCONIn9tFJbgLIyIE80txlR6E0ZFFUZX8zG/ZNbuLUmA20lAqhDjyFsbh2V9WQv4CkhYzKrAuK0e44lyTG7kYjNi0a41QUaNYxY1tkmUxSfg2piGRdb7DrOZgObxUCHqdv+rjJ19JAFLkqTfB/B/Afg9IcQf/gDXfQnAbwM4IoS4oXx2DsDrAD4uhPiu4dxrAPKFEPVPbOwfv9cPD1hXrhChTp9+7JDXyxJKdXVkTFp+k2qTk/QF5uRws8UkZOEi+6/+iqvJn/3ZxzyJybv2/uRPajUEjYe//nWShBMn+JfUh6Eh+hOXlxkwOXs2CTxVVfrqKheaL764LbdrYIDMS9036uxZTSSiCgPfeYdzdGYmwae9XemDesLFi3oZicOHSXPsdsgywePyZRIVs5ks5NlnFXxLJHjC1as6k6mvJ0KVlyOeoEz8+nVdpOB2s/nGRsAsx4jqt2/rJ+TmMui0axeiZjsePaIrUFXcqRWJ2tsBZ1qE9793jz5KSdLUetizB2FXCR71SLh/n0RRCeWgvp6vqbp4kzGi7m6inKqeqawEWluxUd6I3tFUPHqkKwXVPKrmZgJRtq+PbmY1+Qng4DQ1YdXdhP7FQvT1S5ie1l95ejrdoY3Fyyjf6Id5eECX56suzdpaLBU1YDBWjYHxFEwbiI/Vyi425C+gVgwh3TeEpBMkCcJdhoX8BgyjFkOhHEx7pSQwKioUaMwLoVYaQcHaKEzeKa4S1K9WjgtzrgaMSLUYCxdjatZmPIzMDIFGpx81lgmURsZgD03r8TIhIKfYEXA2YMxUg4lYKaaXMxGNSRoY2ywy6jNnUW2egDs+AefaNEwJTduFBMyYTa/FpKUak7ESTG/lIirrfnpJTqAydRbVlkm45SnkbU0jVTK4G9X34HaTnbndRPCPWPWPjxxgSZL06wC+DILMH207HBFCdEuS9AyA3wLwbVA96ADwIoCfB/AtIcTLhvYkAJcAtEJPHP5ZMOH4RSHE937Ifv7wgCUEJ+3LlxmfMiQNq9bVRVefYS/FZPN6udVwaioplVIlQrWlJYaifD5iwtGj27774bC+I6LHQ1plCGYlEsTVS5c4p54/b4g5qc+g1kba2OCsfOZMEoCOjVE4EQrxt3j2LCdfrY1gkMxreFhX8B07lpSgfPky53dZ5sL05EkOl7alyc2b/AuHCZwHD5LiKBXiHzwgPqk5To2NjEOVlirPMDhIeqfGhFwu9mP3biAlBV4vm1e39bBaSbL27VN0KaEQ3YDd3VwEqDlWHR0EgbAN3Yr+QSUb6el6Wb28jDABRAUhQE+S3bUL8ao6DI1Z8OABklxreXnEuZYmGa7VCQJhfz+0on4pKUBDA0RjE2btVegZoKjEGEYqK6PHtDF/Hlm+frJRNb4FkM3X12O1pAGDW+UYGDZjYuJxIKov30KdaQTpM4PspDoOAFBQgFhlHcZt9RheK8LQqFlTOgIkjdWVCTRleOGJjSDNNwKNtgEEs5JSBJwNGBHVGF3Ohddn0sZBCCDHKVCfv4ha0yhKImNI9U/qQTwlD201vxrjtnpMxN0YX3ZiZdWU9Hsodm6hwT4Jj5hA4dYEbMtJqaCIpmbC66jHpKkSE5FizK5nIiFL2i2sUhw1jllUWaZQmphC7uY0rLIOSLKQMGcrxZStBtOiFNPhfKwm0nUGJws45QXUpHpRJnlRFJuGMzHPfd5UkyT+zktL9T+X68cK1D6KgHUFTy65BABTQgiPJEnVAP5fAG0AcgHIYOzqawD+PyFEwniRJEmZYGmmnwRLM/WDpZm+8wH6+cFdguEw41PhMP2ASeoH2o0bdLcdOkRl32OMa36eAoulJWbuGsovAZzgLlygPqCkhN7Ax1JSJiZIq1ZWtF17jf73zU2K8tStO7R9pIzg9fAhb7S2pidGGZSGCwtsY3CQ4HvwoKYup8XjpCZXruh1j06c4I2Uh56aIoCqFYZqa4lv5eVKG+vr9NHducPYWWoqgePgQSA7G7JMgnf9uk4wCgqIT21tyiOHQmRQ3d10Ryl19bB3L1BVhWiMyrzOTn1uT09nN9vbGWPB1BQHq69PjwuVlBDxW1qwFrfjwQOCaSi0La7UCtRUJmAeH+EJasFZgGPS3Aw0N2POXISeXjJCtTQfwLVHUxPQUBFG5uwg+zA6qif1KrI+ua4BM6nV6BuxYWBAj28BZMT19UB90QpK1odgGhogs08kdNrn8SBWVY9xax2G57IxPIwkIHI4GKtqcM2hMjaElOkRLrJkWT+puBgbxTUYt9ZhZL0IY5NmrK/r/bDZgIqyBOozfKiQR+GcH4Xkn9WZHQA4nVjMq8OEpQbj0VKM++1aCBHg6yvNi6A+bZquvo1JWOf9uhwfgEi1I5hdh0lzFSZipZhazcZWOHmFWJS+hjr7NMrFJAojU7CvzXGfNvXrK1nhS6vBtLUKU/ESTG+6EE7YdNyFQK5pETUp0yiTvCiITCM7Pg+TpLexFbdixlYJr9kDr1wCXyQXEXOa1k9JTsAVD6LC5kOZaQYFsRlkJxaYBA8k+2ZLSuiCLCnh7/BHRO34kQOsHxf7UGNYwSD9gC4XtepPyNa/cYOeuIYG5kQ99v0Lh/XK7h0dRLdtQd+ZGbK2xUV6wU6f3qYuFIKT9YUL/EEfO8Y/g/59fZ3A091NPDhxQtlHygik4+NkXn4/Z/OTJ9kn5aR4nLdRscnloutv165tScrvvUcglGWuLI8e1bZvF0qI6PJl3Tvl8fC5NBYXDhM4bt7UqUVtLRFKUREGg8S3Bw+ILSo+7dtH3JXkBEHj7l09FmW1Ejja24HycqxvSHjwgOEm1dOYmqqFm1BaIjjRdneTCW1t6RnBLS1AWxtiecUYHCIYGuNKOTkaRqEwdRlSXy/bMDKh/HygsRFyfSOmwgXoH+BGj0YAKS3lc9WXbyFvaRjS4ACfR40RmUxAWRlEbR3mc+sxOJeDoWFJG1tJ4uurqABqKuKoS52Cc26I7NiImDYbtyMpqcMoqjEyl4XRUb1ChhqKq6wQqM/yo0IeQ0ZwlC5Cow8wMxNRdxUmU+ownijHeNChaTNUy8kBavOWUG2eQGlsAmnBCf2hlUk+kVsAf0YtJqRKTISL4J23G8NiVN7nbKI2lWBUEJlC6lIgCViFyYz5zEpMW6swLdyY2szD0maK9pMBgFRzDDVpPlRYvChOTCN30wtbfBNGW0UmZuw18Jo9mIkXYjbsRFzSQQ2yjDwxh6qUGZTCi4KYD9nxeVgt+kOHY2bMohgz1gr4UAJfNBdrUhaERVE7QiA1vAy3yYcyyyyKZB9cUT8yUyL6z1hVebpcBLOiIoJbYeE/e1xtB7Cekn3olS4ATiB/+Zf8An3uc5przGi9vSxSkZXFGNU2TyDt7l1u0JiSQnVfUnkM/siuX6cGw2wmodL2pFItHqcv7coV/v/hw0QnAwiGw8SUW7f4+961i4q9JLXv2hpPUv16NTUUTWg1nMi+Ll/WsUnZZDfZfbhdRu5yEZ127dJQd0opHqEq8TIzCTwdHYZq8aOjBDAVfOx2TVKO/HwkEmSBnZ1JhcjR1KThE3c07usjGKquRIMbD3V12IpZ0Nub7OkDCKqtrXRNpiXW+EIfPtQrUACcPJqagKYmLIJ5T729Okapc01jI9BQL1BsmSMA9fXpwTNAc+mJ2jr4rB4MjlowOKjXIwQMObTVMmpSppEyOcQBVHOpAH6PqqoQ91RjwlaHkUAGRkb43tT+WCwcm2p3BDWWCeQujUAaH0sGM5MJcLuxWVyNcXMNxjYKMD5pTtJlSBK//tX5q6gyT6BkawzWGQWIVOQEBRhLuTV00cVKMRm0JwG0EIAzW6A2ew6V5ikUx6eQuTAJaX0NRpMdGQhk1mLaUompeAmmVrKxsZXMrtJscdSmz8Jj9qI44aWgIrKRdE7UlAqfoxYzFg+88SJ4N13YkreBmthCdeoMyi0+FCVm4Ar7YBebSZ6K9XgqfLYK+Czl8MlF8EVysWVyACa6MIUsYI+uoMwyizKzD4XCj9zoLDIsWzBbuCIQAlgMp8FvKoHP5IYfhfBHc7FlyYCUovx+ZRmOrXkUwo8y8ywKRACuWAAZ1jBS7ZKebgLwx1NYSHeE+m9u7hPcPf+47QDWU7IPAljXrtHFrm0vv93Gx1mR1ulk5QrDrsOqzc2xuMXiIl10hw8/oa21NWrPh4YYbHjppcfaikSSk3tPn+bcndRWIkGK9957XJE3NRHlDG0JQZbyzjucW7KySKx27zawJhUw3n2XwTWTiRP8sWNJbc3MECvVqkUFBYw7tbQYfiPz8wQeNVvXZiOi7N+v+TzXnqA093gIYg0NCnHc2mLHOzv1hCuHg+673bs1NWF/PxnU5CRvL0mc6NvaCKxWU4LUaLsbLydHCTa1QM7Nx9QUCZKafKsStro6DmttjYBtfpYn9PXpgS+AM3lDA9DQgJCUj4FBCQMDHErV0tJIIuvqgOrcZeYwDQ4SfY3Zs243UFuL9eJajK7mY3hE0piQ+nxZWcyfrXGHUSHGkTozymc0IoOCVLHyakzZajC2kovRMSkJFAHG3CrKEqi1e+GOTyB1dpwvWmUyyuJBLvMgkFGDCeHB5GoOprwmLRyl9quoUKAmO4QK0xSKIpNIDUwm5WkBgEjPwLKrClOWKkzLJZhadWJ+IRmILBagPGcNVdZpuOFFfngaqcuB5D4B2MwowExqNbxSGWZiBZhZy0I0btJOkSTAYYmg2u5DucWHwsQMcje9SIkng1pMNsOf4oHP6oFPFMMXycVSIhNC0vtlFnHkiTl4rD6UmPzIj88iOzqHFKusgYgMExYi6fDZPPCbSjGbKIA/koOYzQGYzewTBOxbiygx+eG2+FEo/HBF/cgwbyIlhccFJKxGUuCPuhCwuhFAIfzxPKwgC3A4IExmSBJgi64jfSOIEpMfRaYg8uQAcs8fgPNkO35Q2wGsp2QfBLCMBWo/8Qmq2J5os7OMT0UizPA1VK5QLR5nW9evP1ZjNtnGxrgFyOIiJ+Lz5w2bYdG2tljItrOTP+ajR9m3x9yG/f10Gy4uJudJGZQhy8uMN3V36zqEJMEEwAn04UOyuKUlXXau6cppgQDzn9Q9CtPSCKr79hnYXCRCRLlzh0t/Fen27mXfUlMhlPBSZydjWWpYp6qK+KSB2Noagae7W/fxbdu7KSFMGB1l99VkW3WSb24mRhUVAdLiAhGqp0dvCyCNbGoCGhoQSXdhZIT4NDzMNYHKWqqqlHJ4tQLpa352XKVJ6m9bpUl1ddgorMKINxVDQ1wbxGL6abm5PK2mMoFykxeWiRHe0NiWCmZVVVjOrcbYVjFGxkyYmNAragD0HFVUAJXuGKqtU8hZHiejUv126onZ2RAVlZjPqsK47MHEfAYmJ5NyeTUtQUXBJqqs0yiNTbDWYiCQ7AO0WiEXl8KfXoMpqRxTWwWY9uulnFRLTwcqXSuosnlRkphGzoYX5lCymw+ShJgzH7P2KkxL5ZiJF8K7mon1jW2gZhbwZC6i0upFCXzIi8zAsRaAJOSk87YsGZi1V8FnLoNPLoRvy4X1WAoE9JWfJGQUW0Mot86iVPIhLz6L7HAANpNBQSJJWIvbMWsth9/ihl8UYjaSg1WRmeSaNyViyIkFUZYSQInkR17CD2ckCEdKHGYT2xGQsLSViiAKELCUwo9CBOMurIgsiFS79ns1x8LI2JpDiTmAYsmPfBGEMzaHdFtUY1uykLCwaUcwko2gpQSVP9GMiuMe/KC2A1hPyT4Ml6AQjEtdv84537D1VLKtr+s79h4+TEqVVFeJZqwxu3cvT3ssHKYKJF5/nbSjtZXyvW1Z+7EYmeC1a/zvhgb2z7ANF215mcir+vNqakgdS0uTTpucJEFTa7ZWVpI1VVdvA7BHjzgg6iRaXk7RRGOj9gPb2KCHUa2tJ4Sy0+0+goUGsHNzdI8+eKDPtqr4obVVk8JPTNC7NzioK6ULC4lNLS2K8DESIVg8eABNKqeWaFIRyuXC8jKBtbc32cuXl6e48BrIEKRZH4G/v58AqwZ48vK0gq2xQjfGJ00YUrx0arFYIQg+tbUcbk/uug4+qkpPtZwcoLoaoroG8+kejEynYGSEwK0+qyRRz1FVxdp9lbYZ2H2jfFk+XzJoZGUBFRUIF1di0lSJsbkMTEzo4hGjN9LjASqdSyiXJ3grD+sAACAASURBVJC9MsU9r7b/XjIzIZeWIZBejSlRhqlVJ6a8Jq3AreE0eIqjqEzxwS1PwbU+BZPfpz+reuP0dKzlVsBrqYBXLoF30wV/yKIp2LXTHAKVmfPwWGbo5gv7qAw0iDEgBGL2TATSKumiSxTCF3ZhYT0F26fVbNsmKu2zcJtmUZCYRc6WD/bYqv7dlljJY0nKgT/FA7+pBLPxfPgjOdhEmnYcAEwigVwRQrnNTwBJ+JEdDcJhiVKkoTQajpkRTOQiYCuDH0UIJnIRimVrTAugCzE1tgZXPAi3xY8CaQ65iSCyEgtIt8ucRoTQASnuwpy1BHNSAYJxF5bkLK4QlVCAORGFfT2E534iDbuOPy4S+8dsB7Cekn2YMSwhyEYuXOAk+fLL3yf2KQRdc2+8wWXkpz9tKLCXfNq9e2RekQgJy3PPPSHxXgiu/N98k6IEt5sop8nu9NOGhti/+XmC4LPPMu6V1Kbq8rt0SQ/cNDXxZEPMSgjO+e+/rxeezc2leGPXLkObQjAYf+sWJ3ZFoqzFnJQAnhCcV+/cIUtRJ2KPh9jU2Ki0KQRR5P59aHtlCMHZta2Nfwo9DQSIwY8e6SUArVaCRFMT/7VZBUGxVxFCqPEaNRCjIlReHkIh9q2/PzkW5XDoBcUrKwHb6jyBcWhIV5MAnIA8HoUi1WAeuRgZZeLv5GRSShLy8xXwqQI8mYuwTY/yvUxM6EVnJYkvsqICqKzEkrMSY0s5mJiUMD6OJMCwWPj6VAAqjU8yrjQxgccCRzk5gMeD5ZxKAtByFianJC3epVpqKr9uHucKyk3TKAxPweqfTq56AXCRUlCAVVcFpk0eeBPF8C6lIxBMztMCFJDMWYXH6kOxPIPc8AwsQZ+u1FTNZsN6Thl8Vg9mUApfNA+zKw5shaXHbl2UsY6KlFmUwEf33OYsq5psW1mu23Iwa/PAby6FXy7AbDgHK1G7BkTqesSBDVSk+lFiDqBA9iMn6kdmfInydcOJG/EUBC0lCFhKEVCYUSiahYQ1Nenelugm8kUQbmsAhVKQYBSZQ3pqnDEtxTZjVgTDWQQjUyGCch7motmIWh18GUo/rbFNODZDKDTN0f2HEJzxENJNm0hL09fJ8tnzMO3rwA9qO4D1lOyfRXQBTpDf+Abnkp/+aebKPNFWVnji2BgnnU996olbKMgyCculS/zddnRQGLHNG0jzegmG09P8Eh85Qka3TZK4qezEcecO57+cHG1LqGTiJ8ucoa9d0wMtFRVsM0lRkawmV9V6zc3EpfJyw6mRCJHk7l3dZeRwEGza2zVXour+u38/WV2el6czJy2LYHmZA/8weX8mlJezE42NQHo6YjGSmL6+5EpGxi2ZamsBe6ogKhkZlNqpzEw9yFRRgbWIDcPDbG9sLHluLSzU67CWFcdh9k4S4UdGkkURZrPmxhOVVQhZizE2QZfl1FRym6orr6IC8ORvIn9jAtLkBNmUEVVU8CkvR7zUg9mUCkwuEXympx/HgMJCoLxMoCJ7CW55CukLU/weqVU6VLNagZISbOWXU74dp3LPZ8AVY6zfXZyAJ20ObniRF/bCFvI93qYkQeS4sJpdhhmLBz5RhJnNHPjnLI8RMIsFKMqJEISkWeTHZ5G17oN5dSmpPQgB2ZaK+YwKzJrdmBVF8MdyEdhIRzSmuw3VbuSlrsFjm0WJOYB8OQBn2I+06DKl7wYUjAgb5mylCFhK6aqL5yIYzkLUpAORJNF9mBFbRFlKUIkbBZETCyIjvgSbTTlVOX89YkVQ5CNkLUYQBQgmcjEfz0LUmp7kzzclYkjbnEehaQ7FZh2MMuVlpNkFf7uSBFlIWN5KQSicgZClCCFTAULChZAKcGlpOP+iOXnLon+i7QDWU7IPBFjXr3PiPXHi+yb9LS7SCzg9TR3BCy/8A6XJJiaYrLWwwMnwxRefKNQQAknbwrvd1E5sI1S0SIT9vH6dM0luLqV7LS2PZTEvGLaESiS27ftkZIpCkA7cuMEZGqC7ob2dyGTIQ4vHOdd3durVgdSiEO3t9DhqQ7exQcDp6tIlcGqJiNZWzvqKi2RujrjU06N7p9Rq5c3NvMRuB8F2eprsqa9PV21IEgdOET/A6UQ0SiAbHOS/xpzV4mJiU00NWYq0tgoNocbHkeSryssjOlVXQ5SVI7hoxcgIyZGq+lZX6jk5ZGSVlUCFOw7HopftjY0RLI1xpIwMDaW2CjyYWHZickrCxAQeE0fY7fw+lJcJeLKXUbA1SaCcmtLjjGq7KSlAWRnkEjcCqR5MxYrhDVjh9SYTL4Df3ZISwJ0fgcc2i6K4F2kLXjLx7UEoiwUoKsKas0xRyxViZiUTvlkpydsJ8LUW5At4MhZQKvlQmPAhe8vPmJUxgKe8u7gzD3NpHvgkMiF/JAfBRWuSzkL9XhVkbKI8JYBiUwD5CT+c0QBS1+YZvzK4DQUkrKbmw28tQ9BcjKCch0DEicWIIymGBdDdl29egNsWRKFEIMqOziE9saILipQ+R+MmzCEfQVspQlIB5gRZ0ToIRMZHs0bW4ZJDKLaGUCDNwSWH2K5JEVno3cVyOBWhRI4ORokczMezEDan8werKBKleAz2rUVkxeZRZJ1HHkJwiXlkxuaR8ann4Tj2gyPWDmA9JftAgCXLVDdcusRJ79/8mydK2AEdZF5/nf99+jRdZ983uX1oiOUlFhc5q7/wwvelaV4vFYLqVlH799N790T2FQoxCNXby47k5/PkJOkebWWFXrx79+h1U/d9OnDgsdCWHpC6d09HkJwcotLu3Y/VPOzpIS4Zi0LU1hKX6uoMC8pYjGPx6JG+PQdANFUk4yrqJRKc63t6eIlxPyaPR9/MNjsbbMfrpdtOzbhVLSNDDyxVVkLYUjA7y9uPjCSHgwzaBlRXA8VFgns2jSpKPCM1Un2HlZVARQWEpwKLghXOJxRyZKzJp9Zh9XgIPmXZq7AHJviiJyeTpeYA2y4rY9Ha3HJMRQoxNWPG9DRJ7PbJXC22UJa3hTLTDLLXvJC803xA1eWo9jstDSgpQSSvFLOWMnjjRfAtklGtrSVPpCpml+TH4EkNoFj4kLPlg2VuVhfSGM1uRzy/GMHUcvhNJfDH8zC7loG5+eRKGGpXXE4Z5Y55uC1U3+XEg0hdCkAKGxQgSkdkWyoW08sQtJbCLwo1JrS2ZUk6FaDirsi+DLc1gCJTELnyHGNO4cWkxGBVcr6GDARTyjBnLkJQzkMo4UQokoW4OSUplgUhYIttoMQSRKElhAIphJxECFlRuui2h7E34zaEolkI2Uowb8pHSHYhFMvGqshIEllIEmCObsERXkC+aR6F5hBcWIAzMY+M+BLSUhL0EiqJ0ZG4GYtbdoTi2VgwF2DBlIuW8x7UHS/GD2o7gPWU7ENzCQ4MUCmRSNCtt622n9GiUeLczZtc3J4+TZD5vuDl9RLpvF4unY8dI/V5gmAjFiMwXr1KjUdaGoHxwIHvs+t4MEj21dPzuHRvW33DeJw4d/u2LkRwOAzVIbZ7MhcW6Mvr7oZW+iA9najU1qbI7yStbVWtp+7SDvB0FZfKygykcHlZl4yrqCdJBJv6ev55PNyPKcH5XQ0pGcsaqZ49BZvoNV1dJdgMDxNFjJO3qhGvqqLKMNWBmRkSotFRvXYgwH+zszV8gscDZJnX2eaEEjsyyt2FIDtVEEp2l8MfYTxqSvHMqXJ11ZxOvWRdWc468sJe1ufzesnQts/4LhdQWgq5uBRBmxveSD5mAhZ4vY/jn3p6cTHgdm2iVPIhP+aDNThDUNvOqCQJyMmBXFiM+bQy+EQxZmN58C+mIBB43P2o5qGV5GzBbfWjSJFspy4HIC3MJwObguCy04WF9DIETCUIyPkIxFwIrqVhbT05bgXwu5KfsQW3LYgi8xzy5CCc0SAcm3MwxaLYbjJMWLIXI2hhfCgkuxCMZGMx6oCMx3OVrIkwSmwhFFvmyFpkBYSMLEvpkCwkLEfsmLcVI2QuRAh5WJCdCEUyETY7kpgWAFjiYWRE5lFoXUC+eQG5IoTsxAIy44tIs8VhterPKgsJK1s2zEczsWApwIIpDwsiBwuJbKzK6ZBT7JBSdXpmiYdh31pEenQRz3wyD43HC/CD2g5gPSX70GNYam2/R4+4hH35ZWXP9yebUYJuNjPcdOzYP1CBZXOTaHTrFiej4mLqzOvqnoh4m5s89dYtds1qpYvv0KEnVpAiU7p7F0nb4rrdBLHm5sf8mevrdCHeu6eHTtTqEG1tnEyTurW2RnB8+FB3eQF6QbympiSl4+oqgaavT8/tVcNdqsihqsrAyFZX6dMbGkrOWZIk0glF7KCC5cqKLsobHyf7U71lmZn6lkmVlQpJXF4mOo2NJVMigAuI0lICjscDuN1Y3rBq7GlqKhksVXwqK1MYlFsgz7oMaXqKCDs9jcdUDkrsCKWlEKVuLKeXYnoxHdPTxKhQKFn1DZDolpQAJcUC7vQlFMRmCDozM1ywqKWaVEtJAYqKIIqKsehg3Gd2Mxu+WQl+f1LhdW1oXS5uxuh2MGcoLzaL1CU/37GaqGakYWlpkPMLsZhWioCpGAE5H/6tbAQXLNraZjv4uJwyStMWGV8SjAU51oOQVleShA6qyXYHFtLcCFpKEJLyMZfIwVw4C0ubKY8JPQDm4BWnLqLYprKgOWTHQkgPL8CMbcAPICFLWDTnIWQtRshUgAW4EIpnYz6cgZjF/rgbUY4jI7aIItsCCkwh5EoLBKHoPNKkrSQQgiQhFpewFLZjwZyPeXM+FqQ8LMrZWIhlYl04uAJVkFGSGNuyh5eQFSfIucxLyBGLyIwtIF1ehT1VICVV0hhXTDZDPnMOKYd28rB+bOyfS3QBgMvtb36T/pjyciZrFXz/1Uw8zrDQ5cucFMrKyL48nn/gHj4fXXzDw/x/t5s6821iCNUiEbribt3S6849UUqumird6+wkvTLW1Nu9m8i0zQ26tUVMevQIWrVwVYne0kJMStoOTFURquKGbRsEar68ggI9QL0OTSauihzUxy0q0nGptFRhZaoLUA0mGUsiqRep6FRWBlitWFmBBjZjY3pOq+quKyoyuOvKWEkBPp/uspueTpb9qX44xW0HtxtLchamp3nq1JSuwTC6wJLqpOZH4YrMQvLN8HlmZh5LtlUVeSgpgSgqxpKjFDPhXPgCZvh8eCLbMZvpHS4qAopzwigxB5AXm4U1NMuxWlhIloobLlKBxy8K4Y/lIriSikAguVsqnlgsdBcWZ2+ixBJEvgjCFQ/CvhKAFJrTx8uIWBYLEq58LKYWa7GlUNyJuU0HlpZNeNLUaLEAhRkbKLbModCkxIJiIaSHQzBvJScCayzfZMNCaglCliLMS3mYl3MQimZhIZKOmHjcmwEh4BDrKLbRHZcrLcApE4DS48uwmMVjP8FowozFeCYWbEVkQnBhQXZiMZqODTgAmw2SSUq6R0pkFdmJBeRbFuGSFpEjFpApLyEjtog0SwwpKQbPgyQhEjNhccuOJTkLi+Y8LEouLMGJpUQmVhIOyDY7JHsqhMmM8+exI7r4cbIPAlgqazl69AlV2Lfb5CSTtebm+Is9d+6xUkvbbWqK7GtqSs/FPXr0+yQUqzY9TRefWl4iI4NotHfvEzbnejIeSRJdWGoS7hNl9D4fqdWjR0QoNc7R2Ejk83iSBkVVoqvVIbZjkurFM2CSfpHKmIzlipRadxoyKaxMvUQNIRmLMAAcepUxeTwGYUYgoCOT16sH+tVZVqVB5eV0qZmtCAQIaKq7bnMzeT7Pzk7eZSI/V6aIQKVD09OPqxrU/ZZKSoCSEiQKihGMu/D/t3fmQZJkd33//DKz7qq+u2d2dlcrodW9gkUgWULWLQhbKyGtwDhw2AgctsMOIGwHDhHGxgiCw1zmCAgiOAwYYTAGLZKMQCYQayFYWbJ2rdWya13Lao+Znemjurqquo48nv/45avMurp7Znqmp1fvE/E6s7LzeC+z6n3z93u/995T54UnU6MoH6xojYpGQw87dw7OnYk551+k1jqv/cQuXNDvnY32sHie3pCbbiJaO8ulgrYfXWiWuXhRb4nt/GyvB/o1OnMGzq5F3BRscpanWQkvEmw9rRabtagsInqtlRWi1TNsFm/WIIR4jYu9BTa3vanbYCmXYWM55ObiplolZpOl8BLV/S283ea4/zV3vXBxja3CTZn7LV5Sl1mnRJzITJGr+EPOlXc462+y7m2zbLZZDLeo97fHRmzP38NB5LNTODPuiosWaYZ17Zdl+1Hl3X3xQK2gQpN1f4dldlgyTRbCHarRHuWSmapPBpFPc1BlR1bZCTZoygpNs0QzbtCK6iTFMqZYGstbEPWp9Js0oiZrwS4r0mSJXRaSJvVol0oQUf6muyi+2gVdnBquRrCiSPvafvSj+sN6y1u0LWduW5Tl0iXtM/W5z+kX2k7HMUNQLEnC2FTvoHr36ldrxTv3mnt76uL75CezTjm33KJ+wTvumOl7tP2r7r9f3XHWBbS6qs1PL33pjI7HoO7ERx5R5XvssUwp6nXN7AtfqFZMMN7gfSE38IOdngpGQ9/xvOfpcswq6/dVYL6Q9k3KPz8rMGnfJM6dGw11s72defMee2y8TchaM1aXbrstfSTDYSYwX/qSqsakiWLnQLr1VnXXrW+wu+fxxBOM0sWL0+66xcWJAbnXQiqtp/WFwKZJtyDoy4Ed8PSmm9irn+NCf5nzF4TzqVHUbk+71kolFZszZ+DMasRNhS024gsUd1Kxefrp8aEw8vd0fR3OnKFTP8tF7yYuxatc7Na5eEnY3MxedvKHep7emo3VmJvKTY18iy6yONwkaG5mltsslpboLZzRKDjZGFk8m+0yrb3ZogPQqCWcqzQ5E2yzLlusJFsshNvUBtv43fbsg4B+oaGWT3CGHW+N7WSJ7eECO8M6/bgw5eKDNGzd7HFTcZs1v8mq7LAUb9OIm9SHTcrecFx8bFh85LMb1dkp6LWaskIzWaQZNdgNa9pPKwjGH54xlIZt6qEK0KrXZJkmC2aXRrRLNd6jUjIUC2bM4hpGHq1+id1hlV1/laas0JIlds0Cf+utG3zlG6e70RyGE6wT4rhcgt2uRup96lP68v+GN2g70ZSLbZIoUvPm3nv1JLWaNmS94hUHjrhsOwH/5V9qxWv7zb7iFWoZjc1UPHngk09qRj/zGfUR2lrlzjtVkWY2bmndYl19+e5DE31rp8Wz3VZxfuQRVYvJdiUb9ZALwgDNmg1mmKVJt96aBTTcfHMuyDGKVCEefVTT+fPjlWKxqAdbX97NN0OxODK0vvSlLE3GFVQqI02yRpD219reVkF76im9v5cujV/TNorlpoowN51jL6mPadOFC5le2Oh7ULE+e1Zv0dmzGqq92D2PXHxaD7pwYWa/JhsIwcYGnDnDYPms9hvqLXJxU2cavnhxul0K9H6urekz3VgO1fJgk8XBJfyti/rmlDeV8wQBrK0RL6+xU1LB2YxTwWkGbG9Pa76lXNb2qrOVFmf8LdbYYilS0SnsbY8PFWLLaW9zpUq7mlk8O2aZ7Ujdes1ukSiWqe+nMfpnudzjTGFH23+8Jouxtv/Uhk1KUXf6e51u6EcBzcIGTX9NxSd1vzWjOq2wShyUZ75RBsmQymCX9UAFaNnbZSlpUo9b1MMmFXpT7j6MIUo89oZldodVmsEaLW+FXVSAWnGDdlxVwSuVkMDPXlqShOKgTbm/S8O0WPVbfOU3Ppvb3zA5O+zhOME6Ia5FG9ZgoPrzV3+lP8oXvECHbDog9iKj21UV+sQn9ETForYVvfKVMzsU59ndVUPq/vuz9oPVVT38zjunRm0aZ2tLAyE+/elMGaxv8MUvHnW6ncQYreQfflj1yPYHsoERz3ue6tFznzsj4j9JtHK3/rt8EAZoJWsV6bbbxsIcrSbZNqannhrXh0JBRcXOYn7LLTn9Hwwyi+nxx6ctJqv+VpVuuUUfXrHI/j4j15zVpnx/Lev9su1BNp09C4V+Ww+yJtD581lMuDVNRDSjZ8+ORtQ2G2doFde5sFXgwoXMEJr3lbX6tLGhbsgzhR112W1f1Ad06dK4uOUr00JhpFLR0hrbhbNsmjU2h4tc2vbZnDCKJgwAGg39zq0vhSOxWUm2aPQ38ZrbenC+nSqP78PyMv3GugqOrKUutgW2ezW2m95U/608lbJhvdHjjL+t7UnSZCneoTHcptxv4vVybx8TGTfi0S6t0QzWNbGsLreowe6wSjvMRGfytolJqMV7bBTV6lmRJkvSohHvUhs2qUTtbNDbCWLj0RqUaXnLmfXDIi2zQCu2gleabttCgyxK/Ra1qMWKr2lR9lgwLepxS92LfkipaLKAjvR7NogDWoMylbveSOP1Luji1HBNgy7QL/fnP58NbO77Kh6vfe2h+qMMhxoWft99mWtobU3bpL7qqw50I4Ie8sAD2txkx+srFnOjij//ACvQxoNbNcq3op85owc///mqCDMa8bpdtZA+9zld5ps2SqVRH1ie8xytXMd+0MboG/wXv5g1FNm2Mks+gOFZzxpT5OFQ77dtX3riCaYqu1pt3FK6+eacqLZaepA1fc6fZ2oCJhvckFels2eJ/SKbm+O6dPHi9AgQoEZXfraHjQ1Yq+6ruFhluphaM7Zzch4rMBsbsL5OsrpOs7DBZrjEpS2PS5f00K2t8fiPfDEWF/UUa2uwvjgcudHqvU1kK1WonZ3pE1iRLZVgdRWzvEK7sjFqy9mKltjulNjeEZrN+Z4/z1MLcnUxYr2wq+1HyTaLSZN6uEOxvaPqPDHo7dh3pVplv7qmlo6XBhkk2pbUHNbY6/pzrw/gS8J6aY/1gorOErssml3qYZNq1KI83MMnmfZ5pgzigJa/QitYpemtquiwyG5UpxXVaIfTbUyjUxlDKepSj2cIT9KiFu1RjjuUi4ZCMB3IESUe7WGJ1qDMnr/Mnr9EiyX2WKCd1GhFNfZNGYolzUPO1XjXXS7o4lRxVYL1qU9prTw1GN987HTvH/2oViIiaoW86lWjuQgPZ3NTTakHH8zapRYXs/5NB0QiAqMRHR56SJe2HgoC7WL0gheoDh3oVtzczML0nnpqfNK+M2eyiLvbbpvZTtbvZ12RZo3UUK1mHrtbb1U9mBou6tKlzErKBzBYS8X3VQmsGp07p7VyLtIw746b17VoeTnTI7tsNNBJIS9dylxyNrjB9tvKNyBVKiO3nDV/zPoG7agycstdTA2gra3p7lOgxVldVRedddWtLQxZNVsUW5v6TGw6SCEWFvQEqci0iutsJqtshwtst4KRFZUPhJj8XpZKasmtrMBqrc+ap9FrS/E21UETb3dHT5LvCW1PZNd9HxYXSRaXaRXXM3dasshO2KC5X6K5Oz0qxqTXs14zrFZ7rAVNVryWtuvE6lqrDZsUuruzxT5H6JfZK62z66+y62kbTytpjKycvbBMZIK5v0+TGMpxl9WgxWrQYsnbY5EWjaRFPdlTaydUS2uqTSvN1yAO2BuUaHuL7PnLtGSJtrdI29RpxXXacYX9pEySs7jG2iiThOKwQ6HfpmH2WA7aLHptFtijbtrUkjbVuE0p6VEsQrEkFN5xF/IKF3Rxargqwer1VHnuu09fn8+e1eCJdDbdo2CMesTuuy+bi7BYzCbCvfXWI4rYrHH0IAsLf/GL50RKZIShGjSfTSegbefap+0ADbbP7ESf4vECXbqUWUaPPZZ1boJs+AYbdXfrrTOtxG43i3GY7ANrsQELNt1000SzXxSphTIrgCEfYlcsZi64vKlTLOoQOLt6/aef1pQPaIDx55N3x62vZ+IShL3MHWctpkuXxufnyN/DajUze1JxiZbW2E6W2Wr6bG1lltP29uy+UcaoyC8vq9CtrMDqimG12GaVbeqD7cxNt72tBbV9siZDA0Ef+vIyLC3Rr6klscMKO9GCRrDtejSbmSU/j1JJn93yQsxaoZVFr8VN6vEu5d4ustucDtef/CGIYBoLdMur7AZr7MrySGh20/aj1n6B/iDXxjVDNwEWSoNMbKTFAns0zB71uEUl3KM83CNIhplbblI1jSFOhI6p0QpWaftL7HlLtFDBacdVWmGVTlwm9ksq1jkX4+h2J4ZitE9p2GbJb7MctFmQNg00VWNN5ahDMUgoFhkbHNcyCD06YYnOoECbBh1/cSR+bVOna6p0kipv/MY6X/umg9oKZuME64Q4Vpfg009rp96HHtLPS0tqOr3sZZc1ZfVgoKHf99+f9WOyzUl33KH6c4gnULEheLMGb61WsxC822+fM45TRqeTRX0/+uh0NPbGRhZZN+GdGyeO1UeWV6NJs6ZYzIIT7HJpabxx3Wge8lp0/vz4xIW2UrIBC1aLzpzRynt0usFAxcOqkTV1rArkK8o0mGCkRjYtL5OIz+7utC7lm20s+XD01dXptFzsErS2VZGsKm1tqWtu0nLKvwwsL2vh0mVYX2bXW9EOp3sFdnbU+NrZyfrhzataKpWRPrG8ZFgptFVcTJN62KS4v6sna6bicqDPzdfCLi4yqK3Q8jVQoGmWaMV1dqM6u90CrT2h252vmZZyGRbrsebJV5FZRK2ZWrxHNWxR6O0hvf3pk+XVSnS+qf3iEu3CirrUvCXaNNhLUndeXKU9LNGLs3akyehLezovDmmgQrPs743EpmY61JI2lahNOe5QSvoUSxMBILk8JUboDAp0wyJtf4m2v0hXGnSkkQlOXKEdVYi8IqZQxASFsRdlET1nEPUphl0Kwy6VuMNS0KHhdanT4TlvfQnPesNz5z+3OTjBOiGuaRtWs6mm0/33Zw0otoPTaK6Mo2Hne7JTW9g6XkRFwrrxch6vg7ENTJ//vCqQHV4AtHK55ZYs2OHWWw/Mq/UQ2r6yjz8+HTxWKIx1L+LcOa0MZ+Z1MFD1yYfONZvTJo2N0c4r0ZkzY3m1wpbXoqefbW8PJgAAHC9JREFUnj6dXa9UxnXIGjgjvQxDRuaNtZKskEz2c7LUauNqlPrSzPIKnbA0MnDsaba3NX/5QXLz5QGt+5eXJ1IjYtlrqeXUambK1GyqOs1qRBtFDshIUFhawiws0isv02SZlllgN1YXXaulp9rdzfR8ntXieWqULSzAQj1mJWizLLssmBYLpkUt3qMybOG1W/qQ8tbmZJikxfcx9QaDylLmNpMF9owKTDuqsBdW2OsVRtOM5Is5K68FP2G11GHJT11oVmSSdiYyUYdi2FU9mGVa5wo+iHw61OkES7T9JbpenQ4qNJ1ELZtOVGY/LpL4RYwfzBRCY8AjIQh7lIZtanRZ9Ds0pENdutToUjMdKkmXctylFHUpBgmFgk6Z4/nThR5GHt1hgW5UojsssPbO17L6xq/icnGCdUJc66CLMWwHpwceUNWxv/ilJTWdXvKSdCjwoyiOkiRqqNiBw/Ph5rY7ku0oe+7cET2VUTpigw12eOKJaWtjZSUbfuHmm1U0ZoxtaBkOVTCefFK16Pz56WH0RDLv3GQgwky97Pczy8gGJ+QtI0s+vHsUXbCu4rG2NhZ9uL+fudrscmsrs0ZmUatl7Tn5tLwM1YrRt/zt7fFkxWRW2Fven7e0lJk5aTILi3SCJXaGdXZbMjJyWi3GXHKz8msfX72u2rSwoMvFRsJS0GHR7Oo0FWEqJK1cmhzIcFZ8eKk0Uqm42qBT0ACAllmgbersxTVaYZV216PdVtfqpDDnAybtaX1f89yoRCwH6i5roO0zNdOhaoUl7ODvt2f3JcufOF+JF2p0giVNXoO2adD1GnSSKu24Sjcu041LdIYlYtSVN8vCyuc3SIYja2bB79LwujSkQ9V0qRq1dEYik/QpFkzmJZwR2GEMDJKCik1YZN/XPHa9BvtSo0uNblKhR5luVKYblYi9gr4pFgpjgRajUxvDW9+S8PJXTo+TeBhOsE6IqxGshx7SCuLOOw8IUDgKdnrbhx/OhgIX0YrUuuxuv/2IfsCMMJzujpSvyIpF1Ztcf9d5g82PY4wW/Mkns0i6p5/O5qq3+P5oJIUxBTrgIv1+pj35QIRZ/YSsB2xSf9bW5oTQ7+xkCmR9dVtb09GHeWy7Tc7NNjJpKjpe3P7+uGWUd7tNxhxMUqtlWrS4ODJwdL0aUhm2tD1ndzdTJGvm5AfdO8ivNzJzNJm6tmm0jFomrYF2xN3bU7FrtbJTz6rzbXlEtC5sNLK0UOxnVorZG1kopcEe0u2oQnU641Els3x/+QvUasTVBt1gUd1jNNSCMTXacYVuXKETlWkPinT3ZWxw/FkW1ZgQeobFUp9FPxUW6VCnQ40ulaRLNelQTvZVWMIuBTMcudrGbsqMCw1jn25cphc06PoLKi5UVVyo0jMaRNGNSuzHJfqJWlsUCqOAisl7b4yG0ReiHoVwn2LYpe73qHv71KVLVXpU2adi9imbHqV4X4Ms4t5ItwqBISik7sgrDBN0gnVCXI1gNZsarGcHIwetgO64Q9PY6OJXQreb9Zr9wheyuZxAf8i28+uzn61WzqG9lMcZDrM+RXY0BuuVGfnkc32K8lFyR2qSiyIVBuuLs2ky0MBebHFx3A9nU70+U1CSRJ/BpP5YDZrnqqrVZrjTljOhGD0zY7SCzatRPuX9snb//LqNLsinkUmziKk36A4LIwPGapJdtlqzyzF5Od/X0+aFo9GARt2wUOzToE092aPYayGdNiNl2tvT8k0OpTSJvbDv682r10dpWGrQlgUNCUhqtJMqnbhKJyzR7mh7VKdzuHDnLwX6wlGrQb0Ushh0WfA6IxeYFZNy1KEUdymFXbxed3ZE4qznkvtfVKiwX1hUMfHqdEXH9Nunyr6p0I3L7Cdl9uMS3bBAPy5iPH+uXk2WgzimEPWoeT0aXpea16Pm7VMzKi4VepQTFZdivE8p7lFI+gS+oeDnrK5ZDx79DfSigqahT48KPa+mSar6mQp9yvRMmV6ibXH9pEjiBdz1zhIvf91R3lLHcYJ1Qhy3S7DbzeYKtBMWWlfX7bdnwwwd2In3KFi1eeyxbLigyZb9fBjduXOqNjPnGJlPkmSak4+Sm+wsC5nXbX097U+UBbjNH33eYhubbNuQVZ/NzYOjxYrFg9VnwpdoDCOLKK89tl2m1Zrf3cdSq41rz4QBQ6ORekcHA8bUyJovrZYKRauVWaUHVbIjf1gjE4vceliqq2DEVdq9YORqs6mTGjaTw1AdZD2Vy1rOsVSOWfBVNOqmra6tpEtp2MbvdxmpU7c7PkL7rIa4yXJaj0K1SlKu0gsa6vZKBaRrMgHpxlrx7ocFuj2P/f3xCNJ57sXJ8hYCQ6PQp+5nYlJVqaIqPUpJTkyifYrRPkHUo0CkI57P8gnOaN8yRq2tXlRQEfHrmaBIlR5l+pTpW0FJRaUXFxmagrbb+QEEwVQo++RlfRMRhD1NUZ9i0lery+tT9fpUpE/Z9CjLgJLps/HWV3D2DS+a/pIfghOsE+J6tWH1++OdaG2Unf2t2jHsbB/Yq3Ixgn6DW62swcgGL1jrZrJ2mozFXl9POxsdvT0tjjOvW779Z3NztkvPZqPRGPe05Ztt6vVDrNTBYDzsLd+gk48MmHfxIBhXnkZj2lyp18cGMe129fRWf/KuNCsS9t1h1u3LVzgjSyIzWkafR8tSRM108Hs5l1pnxnq3q2p7QOU5Xrv5WQaq1VEy1RqDoKZiYSp0kyrdJBWLQcD+vl4qv5w10+88w8aul8uqUZUKVEoJdV8tkbpoO0+ZHhWTikayTzFW4fAH+zpZY78/HuhymP8yv4/vExZS8fDr7EuNvlRGAtIzFfqmRM+U6MUl+kmBXlSkHxfoRwGJ+GMjT0y2ac3abhJDEA8IolRE/AEV6VPzVEAq9NSNx0BT0qOQDCgkfQpxn0I8IPANvqcuPd9nXDhnqHScCIPIpx8FDOKAfhTQjwsMvAoDr8Kz3v7V3PSGF05n/BCcYJ0Q1zXoYg52DLt8H9j8VBZW1FZXx11zGxuHRqMfjX5fVSXfc3Vzc7wj1iR5lcmrzfKyVnxHFLq8181aOnmLxwYQzKsMQC83qTWTmlOrzcnScJi5xiZVx27Ph2zPEwKboVJpXHEmkxWIWg1TrtAb+iPDxGqO1aDuhNGSH4LxoIAKe7+CQC+XGi6j9VEqRNRE2zuqpqttNfE+wSC1kKwi9XrZ51nDZcwyYSZNG9/Xi6YqZUplBn6Vnl+jR2ppmHKW4qIKRlwcubv6ff2q5t9BDhLIuVaID+VCrBaVP1DrQ6ZFo2hUNIpmoIKRpKIR9QkSHdR29BW4DLE0CGHsqYDEBfpedSQgA0oMPLW2+pQZUmRAkUFSZGCKDJIC/bjAMAkYJgGJFyCBr1aY5yN+lilb5qmva5KocMYD3nxXia99zeV5XcAJ1olxIwjWUUgS1ZG8W+7SpXFvWf53Uatl7rhRp9FV1ZPLiKafjTF64XzbTt66mexIM3msNavyPrRZjTCl0qHCZ4xWYlZn8k0z+XRY3558hRYEczVmzCCxQlAuk4U7D4fTSpNPebPE+rKOWuNaisVMdcbMlPFtoV9m31RUAKjQDYvaztFjLO3v67Lf1+WsPmOzsjSJfbEqlTQL+WWlEFGRnEtKnWDqbjMDjZSL1aIIoj7+sIcMB2pB9/vjmTqqGTfLNxgExIWyioNfYSAV+lJhIGVNlOib0pRIDE2BQRIwiAsM4oBB7BMmAVa1Lkc48+tiErxwQJAM8eMhJQaUvSFlb0A5XS/JUO8RQ01mQMEMNSUDAhMSmKGKUBLi++CJwQ8Ez1OB9j0z3unZZsQFXZwuTotgXS7dbtbf1CZrxUyO82qxFY4NeZ7XPtNo5EZGvxKSJPOrHaQy+V7A+QxPVlDGaCU+qTLzkq3Y5xQiiuZrjN1mK/r9/SyCep6hMSvbdr1czjTHrueTrfRHn4tGBzRNNALMD/uZ8vQn1u1nu97vZy7DwzKYZ3J7EGjGSjpE0Gi9VCIplEaVf59MAPqmxMAUxkRgEPnqngp1UFubhkPGIv0Oy9JR1j1PY5JKfkTFH1KSARUZUJJUEBhQSgWhRCoIDCnEA12akCAe4CfhyDrxkxA/CfE8FYjLvpcT68ZAGHuEiU+Y+Awin5ACQ1EnYeiVCD3NYShFhhQYmiIhBd3PFAiNT2jUAgtNoJ+TgGGsFhiehxFPrTLP4663B7z8VfO7o8zjMMG6/DM6vqyx9fZtt13ecUmihsFkbMDjj49rSb7NAua/1HreLGvFo1ptpAmqZ6Dy7ExHLjPYURkOZytKt5tF9FlTwi4nW+rTAgUiLAKLc0O/0Iq6XIZGGdYnzInJ9RkVO6USiV9gMJQpbRkMdN02z9ltuo8wGBTTtHiktqOp9Xxdap8ZEPiazXmpUEiL4UfqMksrd/vGXzQDbXNJ3/4X4yZBom/+MkyVaDjUwoRh9nlWIbD3Pvd5TATy+x2yXiyQ+AVCv8xQtOIfispTKMW08teKf0iRtmkQSqAVv2jFPzSBbvN8QnyG+PoZFZjYpCIgMp2L/P3Ply537yUtqvFB0treJAbPxPgmQqKQwIR4wyF+ElL0IgqEFCWkIBFFCalKlwIhgYQUvJAA3cc3EUGgsuYTqRUWRvhJyEL/jcDldxw+DCdYKSJSB34U+HvAEvDXwA8ZYz5wohl7hpAfneDWy58mZ4ooyrQib7Hs72cxIHkt6fXGp7uf9IbNEkjdXkSkiOctzdSLchlKjWntsOu2Us6v54Z7G8e6/qzSWDWx6/bz3h5TpkPuszccUgEqV+JLmrwZttozJussalXmiCn2UtdX+tau64XRG/vQFAgTn2ES0I3K7MZVwkhGFpHVoPzSpnlGx9h2Jip0DpaiUfEntotAkPY5CryEkq8V+qiCzy0LEhEYa6NoqrCPbyt7QgKJ8D2t4H3SdUI8ifAlxEsivDjEiyI8EjxR95tIGhQxrxRyUOkAL/dF97P9TWKI8YkSjzASImPXPWLjEcVCSIFYAiKvQESBSAoMqRAREIr9X8CLaLAy4yt+tTjByrgHeBnwbuBvgG8H7hGRtxljPnSSGXNMY4Pw5g6We8zE8bhmTK7bl/pud1pLJl/+5zebCCIloIQx4/0TrFvVHjP7eF36/vhABPOWk+sHJt9kb9ZJmLZvpCkZ4sXpm3ocjitKv48fRVTCkEoYauHt/ybX7efLNe+Ysf1AC4kD/ifZR8/Twvs+BBrBF6WVspWh2CsQERARaGUtBUITaMWfWlGx8eiZMqGpE+MTo9ZUjKeigE8sKgqhXfc8IjTFXiYYceSRGFH3m6dfCiNeakkdKGPZ+gyLzN4SEcBLLbJZtzXWEkgS4yUREkejz/n1DVa4TCfMkXCCBYjIW4A3A+80xtyTbvtz4CuAnwacYH2ZYyO1jzRixwkTx7P1wK7HcfY5rxNxrAI8a584FsKwQBwXRsfk94ljTfnod7jKdS7LQXfkddLPnkcWPOCPf/bEEHgJPrHKkagLLRArOblkooltkbrciPFMTJGIisR4aPJNjCf6Py/dLibRdS/9nGSfJYnxSBAvRkhSsdBjJI4RE+NhMstLMits0iLTezvHIjuKS9RL97fNtPMeXv0u4AomxDoEJ1jK3UALeL/dYIwxIvKbwC+LyIuNMQ+fWO4cjsvAVsCO+dhglrzQ2qWuC0niE8c+cVwc29eYWfuPrxsDcQLDGdsPW7f722vO2s/+f9657OfJ7XA0Yc9WDlg/4H/XRq6cYFnuAB42xkzOY/Bg/v/XN0sOh+NaYV2sNsrPcTpwgqWsAp+bsX0n9/8xROSwePWrHSTJ4XA4HDmuZgjVZxrmCv/ncDgcjuuAs7CUbWZYUTCKzNyZ/Me8jm2W1AJzVpbD4XAcE87CUv4aeJGITN6Pl6bLh65zfhwOh8MxgRMs5R60s/DbJrZ/G/BZFyHocDgcJ49zCSofAv4c+DURWUU7Dr8L+NvA208yYw6Hw+FQnGAx6nP1DnRoph9Fra2H0Y7EHzzRzDkcDocDwI3Wfq0QkQSQxaueAtjhcDi+PGi1WqA2xMzmKidY1wgRidA2wr0rONyqXOv4cnTD48r8zOfLrbzgyny5LACJMWam988J1g2I7ZR8WOj8MwlX5mc+X27lBVfm4z63ixJ0OBwOx6nACZbD4XA4TgVOsBwOh8NxKnCC5XA4HI5TgRMsh8PhcJwKnGA5HA6H41TgBMvhcDgcpwLXD8vhcDgcpwJnYTkcDofjVOAEy+FwOBynAidYDofD4TgVOMG6johIXUR+XkQuiEhPRP6PiHzjEY57j4iYGenp65HvK0VEbhGRnxORj4lIJ83z6y/j+K8RkT8Tka6INEXkd0Xk5muY5avmasosIr8x5zl//Bpn+6oQkTelef+siOyLyJMi8j4ReenhR4OIPFdE/lBEWiLSFpEPiciLr3W+r5SrKe8p/i1/nYh8WESeEpG+iGyKyEdE5O8e8fhjecZuPqzryz3Ay4B3o5NEfjtwj4i8zRjzoSMc//VAJ/d5eOw5PF5uB74VuB/4M+BQcbaIyIuAe4FPAt8M1IAfAe4Vka82xnQOOPwkueIyp3TQ55ynfQz5upb8c2AV+BngEeAM+h3/pIi83hgzV3BFZAP4C+ASOmlqBPx74H+lz/nJa535K+CKy5vjtP2Wl4HPAr8OPJ1+/mfAh0TkW40xvzvvwGN9xsYYl65DAt4CGODu3DYBPgY8csix70mPXTrpclxmmb3c+jvSMrz+iMf+HnAeqOW2vRCIge896bJdozL/BrB70mW4gjJvzNi2BDSBPzjk2J8AesC53LZVdFqeXzrpsl2D8p7K3/KcsgTAE8BHrtczdi7B68fd6Pww77cbjD653wReeCO7QK4UY0xyJceJSAF4K/D7xphu7nz/D/g48E3Hk8Pj50rLfJoxxlyasW0X+DxwyyGH3w38qTHmfO7YbeCDwDuPM5/HxVWW9xmDMSZC67TwkF2P7Rk7wbp+3AE8PKNCezD3/8N4RETitA3sV1JT+5nIVwAV4KEZ/3uQo92r00pdRC6mz/lLIvLTIlI/6UxdLiKyjj6nWc/Q7lMBnjtnnweBjdPyHT9KeSc4lb9lEfFEJBCRcyLyg8DzUdfovP2P9Rm7NqzrxyrwuRnbd3L/n8cXge8DHkB93a9GfeZvEpGvMcY0jzOjNwD2XuzM+N8OUBGRijGmdx3zdD34NPB/0R+3j7ZzfDfwGhF5tTHmsDfZGwIREeCX0Rfinzpg12XULT7vOYN+F6YsmhuJyygvnP7f8u+ReTj2gG8xxvzJAfsf6zN2gnV9OWhYkbn/M8b81sSmj6SRY/8T+E7gh48hbzciV3S/TivGmMk31Q+LyGfRyvDvA++9/rm6In4Sbb/7DmPMI0fY/7Q/5yOX9xnwW3438OPAWeAfAL8nIu8yxvzOIccdyzN2LsHrxzazraiVdDnrDWQuxpg/BS4Ar7rKfN2IbKfLeferZ4zpX8f8nCTvBRJOyXMWkR8Bvgf4l8aY3zhk9yZaWR3b7+J6c5nlnclp+i0bYx41xnzSGPNBY8y3Ah8GflFE5mnJsT5jJ1jXj78GXjTjwdq+G0f1fefx0MrsmcajaFTRrLaql3Jl9+q0Iunyhn/OIvJDqLvr3caYnz9s/9Sl+yjzn/PmrACHG4XLLe8hnNbf8idQt9/6rH8e9zN2gnX9uAcNfX3bxPZvAz5rjHn4ck4mIt+A9v+4oTuVXglpW80fAd8kIlW7XUSej76Fvu+k8nYC/EP0d3pDP2cR+QHg+4HvN8b85GUceg/w9SJyNneuFfR3csM+56so76xzncrfctp293pgl8wrMotje8ZutPbrRPpw/wz4SrKOw+9CBevtxpgPpvvdC7zOGCO5Yx8A/gvacS8Evg74N2gHvpenIbU3JCLyzenqy9Fyvwe1NrvGmD9O93kMwBjz7NxxL0bf3j6ONmTbjsMF4E5jzA3bmfZKyiwitwG/BfwO2jDvA28Gvgv4FPDaNIz4hkNEvgd9Rv8DfUZ5BsaYB9L97mX6u30GDTY5D/wgWafS5wNfbYx5/JoX4DK5yvKeyt+yiPw28CX0u7gF3ITWX38H+G5jzC+k+93LtXzGJ9357MspAQvAL6Bfzj46GsI7Jva5l7SLVm7b76B9PLpoZNEX0VDSlZMu0xHKbOakx3L7PJb/nNv+cuAjabl30QilW0+6TNeizKhb5X3p9l76/Xg4/YFXTrpMh5T33iOWeeq7nW5/Hto/cQ8d/eGPgZecdLmuRXlP628ZfXG6D7WkonT5YeBts+7NtXrGzsJyOBwOx6nAtWE5HA6H41TgBMvhcDgcpwInWA6Hw+E4FTjBcjgcDsepwAmWw+FwOE4FTrAcDofDcSpwguVwOByOU4ETLIfD4XCcCpxgORyOK0ZEflFEnkqHHnM4rilOsBwOxxWRitTbgfcbN2SO4zrgBMvhcFwpLwduBv7wpDPi+PLACZbDcYMgIv9ZRIyITM0tJCK3i8hQRH7pKs7/7en53yQi/0FEviQiPRH53yLyynSf14nIx0SkKyIXROT7Dzjl3UAL+PP02LKIvEdEPisi+yKyKyKfEZGrmn7D4bAEJ50Bh8Mx4sF0eQepCOT4cXQE9x84huv8R3T6kp8DiuiMuR8WkXcBvwb8MvDbwLcAPyQif2OMee+M89wN/JHR+csAfhH4x+j0GT+TXuN5wBuPIc8OhxMsh+MGYqZgicjXAe8Evs8czwy8PvBKY8wwPf/D6NQPvw+8yhjzyXT7r6FzIH0nMCZYIvIi4AXovEaWu4E/Nsa86xjy6HBM4VyCDseNw6fT5eR04j8FPIFaLcfBL1mxSvmLdPlxK1YA6T6fQK2kSd4BDIA/yW1rAS8RkVnToTscV40TLIfjBsEYsw1cAF5it4nItwCvAv6tMaafbgtE5OdEZCdtJ/pVESldxqUenbhuM139mxn7NoHVGdvvBv7UGNPJbftX6ESUnxGRL6b5eruIuHrGcSy4L5LDcWPxIKlgiUgR+DHgk8B/ze3zfcDrUEvseen+P3YZ14gvc/sYInIL8LVMRAcaY94PPBv4R+hM0W9K97k3LYvDcVU4wXI4biweBJZSUfhO4CuA75no5/RPgB82xpw3xmwC7wG+4zpaMu9Ap4P/wOQ/jDE7xpj3GmP+KZr3nwBeg/bXcjiuCidYDseNhQ28eA0a0PA+Y4xtY0JEloBbgQdyx9wP2O3Xg7uBj6ViafPlp3kbkYqszefKdcqb4xmMixJ0OG4srGD9LNAAvnfi/4102cpt25343zVDRJaB1wLvnpGvCyLyAVSkLgHPAf4F2g72wWudN8czHydYDseNxSPAENgAftYY84WJ/7fT5SKwla4vTfzvWvI2tN6YHN1iHxXZNwFvBupoAMkHgB8zxpy/DnlzPMMRNwSYw3G6EJHHgX9tjPmD9PM3AP8NWDXGJNf42vcAzzHG3Hktr+NwzMJZWA7H6eNXgX8nIvcBIRp08evXWqxS7gN+5Tpcx+GYwllYDscpQ0QC4D+h4eMe8N+B77L9tByOZypOsBwOh8NxKnBh7Q6Hw+E4FTjBcjgcDsepwAmWw+FwOE4FTrAcDofDcSpwguVwOByOU4ETLIfD4XCcCpxgORwOh+NU4ATL4XA4HKeC/w85WAKsrXZX1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "v_0 = np.linspace(0.5, 3., 100)\n",
    "for k in range(1,20):\n",
    "    omega_01 = 9.8 * (2 * np.pi * k + 0.5 * np.pi) / (2 * v_0)\n",
    "    omega_02 = 9.8 * (2 * np.pi * k + 1.5 * np.pi) / (2 * v_0)\n",
    "    ax.plot(v_0, omega_01, 'b', lw=0.5)\n",
    "    ax.plot(v_0, omega_02, 'r', lw=0.5)\n",
    "ax.set_xlabel(r'$v_0$ m/s')\n",
    "ax.set_ylabel(r'$\\omega_0$ rad/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we discussed, the uncertainty rises from our inability to perfectly control the initial conditions of the coin toss experiment.\n",
    "Let us suppose that a typical human throws the coin with an initial velocity:\n",
    "$$\n",
    "v_0 = 2.5 \\pm 0.2\\;\\text{m/s},\n",
    "$$\n",
    "and an initial angular velocity of\n",
    "$$\n",
    "\\omega_0 = 400\\pi \\pm 50\\;\\text{rad/s}.\n",
    "$$\n",
    "For the time being let's interpret the $\\pm$ as if it means that all values plus or minus that value are equally likely.\n",
    "\n",
    "Let us now assume that the experiment is repeated 1,000 times and that each time the initial conditions are drawn randomly.\n",
    "Then we are going to count the number of times we get an H.\n",
    "Dividing by the total number of random experiments, we will get the frequency of H's.\n",
    "The code is given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure_freq_of_H(v_0_lower=2.3, v_0_upper=2.7, \n",
    "                      omega_0_lower=200*np.pi-100, \n",
    "                      omega_0_upper=200*np.pi+100,\n",
    "                      N=10000):\n",
    "    num_Hs = 0.\n",
    "    for n in range(N):\n",
    "        v_0 = np.random.rand() * (v_0_upper - v_0_lower) + v_0_lower\n",
    "        omega_0 = np.random.rand() * (omega_0_upper - omega_0_lower) + omega_0_lower\n",
    "        num_Hs += 1 if X(v_0, omega_0) == 'H' else 0\n",
    "    return num_Hs / N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run this code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.503"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "measure_freq_of_H()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is our first Monte Carlo simulation. As we will see in a few lectures, the empirical frequency that we measure this way converges to the probability of the coin turning out heads in the limit of $N\\rightarrow\\infty$.\n",
    "This is known as the *law of large numbers*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Run ``measure_freq_of_H()`` two or three times. What do you observe?\n",
    "+ If our assumptions about the physics of the coin toss are correct, then the frequency we get from the Monte Carlo should be very close to the probability of getting heads. That should be very close to 0.5. The Stanford statistician Persi Diaconis claims that he is able to get any outcome he likes from a coin toss experiment. In what regime of the initial condition space do you think he operates? Remember that: 1) he needs to be able to control for the outcome; but 2) he also needs to full you that the is throwing the coin regularly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The common sense assumptions that give rise to the basic probability rules.\n",
    "\n",
    "> Probability theory is nothing but common sense reduced to calculation. Pierre-Simon Laplace, Théorie analytique des probabilités (1814)\n",
    "\n",
    "Consider the following three ingedients:\n",
    "+ A: a logical sentence\n",
    "+ B: another logical sentence\n",
    "+ I: all the information we know\n",
    "\n",
    "No other restriction apart that A and B are not contradictions.\n",
    "\n",
    "We need a bit of notation so that we write less math:\n",
    "\n",
    "+ $\\text{not}\\;A \\equiv \\neg A$\n",
    "+ $A\\;\\text{and}\\;B \\equiv A, B \\equiv AB$\n",
    "+ $A\\;\\text{or}\\;B \\equiv A+B$\n",
    "\n",
    "Now, let's try to make a robot that can argue under uncertainty.\n",
    "It should be able to take logical sentences (such as $A$ and $B$ above) and argue about them using all the information it has.\n",
    "What sort of system should govern this robot.\n",
    "The following desiderata seem reasonable:\n",
    "\n",
    "+ Degrees of plausibility are represented by real numbers.\n",
    "+ The system should have a qualitative correspondence to common sense.\n",
    "+ The system should be consistence in the sense that:\n",
    "   - If a conclusion can be reached in two ways, each way must lead to the same result.\n",
    "   - All evidence relevant to a question should be taken into account.\n",
    "   - Equivalent states of knowledge must be represented by equivalent plausibility assignments.\n",
    "   \n",
    "[Cox's theorem](www.sciencedirect.com/science/article/pii/S0888613X03000513?via=ihub) shows that:\n",
    "\n",
    "> The desiderata are enough derive the rules of probability theory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Talking about probabilities\n",
    "\n",
    "We read $p(A|BI)$ as:\n",
    "\n",
    "+ the probability of A being true given that we know that B and I are true; or\n",
    "+ the probability of A being true given that we know that B is true; or\n",
    "+ the probability of A given B."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpratation of probabilities\n",
    "\n",
    "The probability $p(A|BI)$ is a number between 0 and 1 quantifying the degree of plausibility that A is true given B and I.\n",
    "Specifically:\n",
    "\n",
    "+ $p(A|B,I) = 1$ when we are certain that A is true if B is true (and I).\n",
    "+ $p(A|B,I) = 0$ when we are certain that A is false if B is true (and I).\n",
    "+ $0< p(A|B,I) < 1$ when we are uncertain about A if B is true (and I).\n",
    "+ $p(A|B,I) = \\frac{1}{2}$ when we are completely ignorant about A if B is true (and I)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The rules of probability\n",
    "\n",
    "There are two rules of probability from which everything else can be derived.\n",
    "These are direct consequencies of the desiderate and Cox's theorem.\n",
    "They are:\n",
    "\n",
    "+ The **obvious rule**:\n",
    "$$\n",
    "p(A|I) + p(\\neg A|I) = 1.\n",
    "$$\n",
    "The sum rule is obvious. It states that either $A$ or its negation $\\neg A$ must be true.\n",
    "(It is vitally important that you do not try to apply probability in a system that includes contradictions.)\n",
    "\n",
    "+ The **product rule** (or Bayes' rule or Bayes' theorem):\n",
    "$$\n",
    "p(A,B|I) = p(A|B,I)p(B|I).\n",
    "$$\n",
    "The product rule is not obvious.\n",
    "Understanding it requires a bit of meditation.\n",
    "It states that the probability of A and B is the probability of A given that B is true times the probability that B is true.\n",
    "Even though the correspondance is not one to one, visualizing events using the Venn diagrams helps in understanding the product rule:\n",
    "![Venn diagram](venn.png)\n",
    "In this diagram:\n",
    "    - $p(A,B|I)$ corresponds to the brown area (normalized by the area of I).\n",
    "    - $p(B|I)$ is the area of $B$ (normalized by the area of I).\n",
    "    - $p(A|BI)$ is the brown area (normalized by the area of B)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Drawing balls from a box without replacement (1/3)\n",
    "\n",
    "Consider the following information I:\n",
    "\n",
    "> We are given a box with 10 balls 6 of which are red and 4 of which are blue. The box is sufficiently mixed so that when we get a ball from it, we don't know which one we pick. When we take a ball out of the box, we do not put it back.\n",
    "\n",
    "![Urn](urn.png)\n",
    "\n",
    "Now, let's draw the first ball.\n",
    "Here is the graphical causal model up to this point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Urn1 Pages: 1 -->\n",
       "<svg width=\"222pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 221.82 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>Urn1</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 217.82,-112 217.82,4 -4,4\"/>\n",
       "<!-- reds -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>reds</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"46.94\" cy=\"-90\" rx=\"46.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"46.94\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\"># red balls</text>\n",
       "</g>\n",
       "<!-- first -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>first</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"104.94\" cy=\"-18\" rx=\"40.13\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">1st draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;first -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>reds&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M60.69,-72.41C67.91,-63.69 76.9,-52.85 84.89,-43.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"87.61,-45.4 91.3,-35.47 82.22,-40.94 87.61,-45.4\"/>\n",
       "</g>\n",
       "<!-- blues -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>blues</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"162.94\" cy=\"-90\" rx=\"50.76\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"162.94\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\"># blue balls</text>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;first -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>blues&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M149.2,-72.41C141.98,-63.69 132.99,-52.85 125,-43.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"127.67,-40.94 118.59,-35.47 122.28,-45.4 127.67,-40.94\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a22f58d90>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gu1 = Digraph('Urn1')\n",
    "gu1.node('reds', label='# red balls', style='filled')\n",
    "gu1.node('blues', label='# blue balls', style='filled')\n",
    "gu1.node('first', label='1st draw')\n",
    "gu1.edge('reds', 'first')\n",
    "gu1.edge('blues', 'first')\n",
    "gu1.render('urn1_graph', format='png')\n",
    "gu1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's say that we draw the first ball.\n",
    "Let B_1 be the sentence:\n",
    "\n",
    "> The first ball we draw is blue.\n",
    "\n",
    "What is the probability of $B_1$?\n",
    "Our intuition tells us to set:\n",
    "$$\n",
    "p(B_1|I) = \\frac{4}{10} = \\frac{2}{5}.\n",
    "$$\n",
    "This is known as the *principle of insufficient reason*.\n",
    "We can now use the **obvious rule** to find the probability of drawing a red ball, i.e., of $\\neg B_1$.\n",
    "Of course, $\\neg B_1$ is just the sentence:\n",
    "\n",
    "> The first ball we draw is red.\n",
    "\n",
    "So, let's call it also $R_1$.\n",
    "It is:\n",
    "$$\n",
    "p(R_1|I) = p(\\neg B_1|I) = 1 - p(B_1|I) = 1 - \\frac{2}{5} = \\frac{3}{5}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the graphical model representation after we observe the first draw?\n",
    "We need to fill the node corresponding to the first draw with color:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Urn3 Pages: 1 -->\n",
       "<svg width=\"227pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 226.82 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\n",
       "<title>Urn3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-184 222.82,-184 222.82,4 -4,4\"/>\n",
       "<!-- reds -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>reds</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"46.94\" cy=\"-162\" rx=\"46.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"46.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\"># red balls</text>\n",
       "</g>\n",
       "<!-- first -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>first</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"104.94\" cy=\"-90\" rx=\"40.13\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">1st draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;first -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>reds&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M60.69,-144.41C67.91,-135.69 76.9,-124.85 84.89,-115.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"87.61,-117.4 91.3,-107.47 82.22,-112.94 87.61,-117.4\"/>\n",
       "</g>\n",
       "<!-- second -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>second</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"104.94\" cy=\"-18\" rx=\"43.02\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">2nd draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;second -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>reds&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.54,-143.77C44.73,-125.21 45.48,-95.22 55.94,-72 61.04,-60.69 69.53,-50.2 78.01,-41.59\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"80.58,-43.98 85.39,-34.53 75.74,-38.92 80.58,-43.98\"/>\n",
       "</g>\n",
       "<!-- blues -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>blues</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"167.94\" cy=\"-162\" rx=\"50.76\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"167.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\"># blue balls</text>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;first -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>blues&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M153.01,-144.41C144.98,-135.48 134.92,-124.31 126.09,-114.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"128.63,-112.08 119.33,-106.99 123.42,-116.76 128.63,-112.08\"/>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;second -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>blues&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M168.15,-143.72C167.73,-125.12 165.13,-95.09 153.94,-72 148.45,-60.66 139.71,-50.05 131.13,-41.35\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"133.35,-38.63 123.71,-34.22 128.5,-43.68 133.35,-38.63\"/>\n",
       "</g>\n",
       "<!-- first&#45;&gt;second -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>first&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M104.94,-71.7C104.94,-63.98 104.94,-54.71 104.94,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"108.44,-46.1 104.94,-36.1 101.44,-46.1 108.44,-46.1\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a23042e50>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gu3 = Digraph('Urn3')\n",
    "gu3.node('reds', label='# red balls', style='filled')\n",
    "gu3.node('blues', label='# blue balls', style='filled')\n",
    "gu3.node('first', label='1st draw', style='filled')\n",
    "gu3.node('second', label='2nd draw')\n",
    "gu3.edge('reds', 'first')\n",
    "gu3.edge('blues', 'first')\n",
    "gu3.edge('first', 'second')\n",
    "gu3.edge('reds', 'second')\n",
    "gu3.edge('blues', 'second')\n",
    "gu3.render('urn3_graph', format='png')\n",
    "gu3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the sentence $R_2$:\n",
    "\n",
    "> The second ball we draw is red.\n",
    "\n",
    "What is the probability of $R_2$ given that $B_1$ is true?\n",
    "We just need to use common sense to find this probability:\n",
    "+ We had 10 balls, 6 red and 4 blue.\n",
    "+ Since $B_1$ is true (the first ball was blue), we now have 6 red and 3 blue balls.\n",
    "+ Therefore, the probability that we draw a red ball next is:\n",
    "$$\n",
    "p(R_2|B_1,I) = \\frac{6}{9} = \\frac{2}{3}.\n",
    "$$\n",
    "\n",
    "Similarly, we can find the probability that we draw a red ball in the second draw given that we drew a red ball in the first draw:\n",
    "+ We had 10 balls, 6 red and 4 blue.\n",
    "+ Since $R_1$ is true (the first ball is red), we now have 5 red and 4 blue balls.\n",
    "+ Therefore, the probability that we draw a red ball next is:\n",
    "$$\n",
    "p(R_2|R_1,I) = \\frac{5}{9}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's consider a second draw without observing the result of the first draw.\n",
    "What is the graphical causal model now?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Urn2 Pages: 1 -->\n",
       "<svg width=\"227pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 226.82 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\n",
       "<title>Urn2</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-184 222.82,-184 222.82,4 -4,4\"/>\n",
       "<!-- reds -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>reds</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"46.94\" cy=\"-162\" rx=\"46.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"46.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\"># red balls</text>\n",
       "</g>\n",
       "<!-- first -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>first</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"104.94\" cy=\"-90\" rx=\"40.13\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">1st draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;first -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>reds&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M60.69,-144.41C67.91,-135.69 76.9,-124.85 84.89,-115.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"87.61,-117.4 91.3,-107.47 82.22,-112.94 87.61,-117.4\"/>\n",
       "</g>\n",
       "<!-- second -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>second</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"104.94\" cy=\"-18\" rx=\"43.02\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">2nd draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;second -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>reds&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.54,-143.77C44.73,-125.21 45.48,-95.22 55.94,-72 61.04,-60.69 69.53,-50.2 78.01,-41.59\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"80.58,-43.98 85.39,-34.53 75.74,-38.92 80.58,-43.98\"/>\n",
       "</g>\n",
       "<!-- blues -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>blues</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"167.94\" cy=\"-162\" rx=\"50.76\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"167.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\"># blue balls</text>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;first -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>blues&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M153.01,-144.41C144.98,-135.48 134.92,-124.31 126.09,-114.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"128.63,-112.08 119.33,-106.99 123.42,-116.76 128.63,-112.08\"/>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;second -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>blues&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M168.15,-143.72C167.73,-125.12 165.13,-95.09 153.94,-72 148.45,-60.66 139.71,-50.05 131.13,-41.35\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"133.35,-38.63 123.71,-34.22 128.5,-43.68 133.35,-38.63\"/>\n",
       "</g>\n",
       "<!-- first&#45;&gt;second -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>first&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M104.94,-71.7C104.94,-63.98 104.94,-54.71 104.94,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"108.44,-46.1 104.94,-36.1 101.44,-46.1 108.44,-46.1\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a22feae90>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gu2 = Digraph('Urn2')\n",
    "gu2.node('reds', label='# red balls', style='filled')\n",
    "gu2.node('blues', label='# blue balls', style='filled')\n",
    "gu2.node('first', label='1st draw')\n",
    "gu2.node('second', label='2nd draw')\n",
    "gu2.edge('reds', 'first')\n",
    "gu2.edge('blues', 'first')\n",
    "gu2.edge('first', 'second')\n",
    "gu2.edge('reds', 'second')\n",
    "gu2.edge('blues', 'second')\n",
    "gu2.render('urn2_graph', format='png')\n",
    "gu2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find the probability that we draw a blue ball in the first draw (A) and a red ball in the second draw (B).\n",
    "We have to use the **product rule**:\n",
    "$$\n",
    "p(B_1, R_2|I) = p(R_2|B_1,I) p(B_1|I) = \\frac{2}{3}\\frac{2}{5} = \\frac{4}{15}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other rules of probability theory\n",
    "\n",
    "All other rules of probability theory can be derived from the two basic rules.\n",
    "Here are some examples.\n",
    "\n",
    "#### Extention of the obvious rule\n",
    "For any two logical sentences $A$ and $B$ we have:\n",
    "$$\n",
    "p(A + B|I) = p(A|I) + p(B|I) - p(AB|I).\n",
    "$$\n",
    "In words: the probability of A or B is the probability that A is true plus that probability that B is true minus the probability that both A and B are true.\n",
    "This is very easy to understand intuitively by looking at the Venn diagram:\n",
    "![Venn diagram](venn.png)\n",
    "The probability $p(A+B|I)$ is the area of the uninion of A with B (normalized by I).\n",
    "This area is indeed the area of A (normalized by I) plus the area of B (normalized by I) minus the area of A and B (normalized by I) which was doublecounted.\n",
    "\n",
    "Let's see a formal proof of this.\n",
    "$$\n",
    "\\begin{split}\n",
    "p(A+B|I) &=& 1 - p(\\neg (A+B)|I)\\\\\n",
    "&=& 1 - p(\\neg A, \\neg B|I)\\;\\text{(obvious rule)}\\\\\n",
    "&=& 1 - p(\\neg A|\\neg B, I)p(\\neg B|I)\\;\\text{(product rule)}\\\\\n",
    "&=& 1 - \\left[1 - p(A|\\neg B, I)\\right]p(\\neg B|I)\\;\\text{(obvious rule)}\\\\\n",
    "&=& 1 - p(\\neg B|I) + p(A|\\neg B, I)p(\\neg B|I)\\\\\n",
    "&=& 1 - p(\\neg B|I) + p(A\\neg B|I)\\;\\text{(product rule)}\\\\\n",
    "&=& 1 - p(\\neg B|I) + p(\\neg B|A,I) p(A|I)\\;\\text{(product rule)}\\\\\n",
    "&=& 1 - p(\\neg B|I) + \\left[1 - p(B|A,I)\\right]p(A|I)\\;\\text{(obvious rule)}\\\\\n",
    "&=& 1 - p(\\neg B|I) + p(A|I) - p(B|A,I)p(A|I)\\\\\n",
    "&=& 1 - \\left[1 - p(B|I)\\right] + p(A|I) - p(B|A,I)p(A|I)\\;\\text{obvious rule})\\\\\n",
    "&=& p(A|I) + p(B|I) - p(B|A,I)p(A|I)\\\\\n",
    "&=& p(A|I) + p(B|I) - p(AB|I)\\;\\text{(product rule)}.\n",
    "\\end{split}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The sum rule\n",
    "\n",
    "This is the final rule we are going to consider in this lecture.\n",
    "It is one of the most important rules.\n",
    "**You absolutely have to memorize it.**\n",
    "It goes as follows.\n",
    "\n",
    "Consider the sequence of logical sentences $B_1,\\dots,B_n$ such that:\n",
    "+ One of them is definitely true:\n",
    "$$\n",
    "p(B_1 + \\dots + B_n|I) = 1.\n",
    "$$\n",
    "+ They are mutually exclusive:\n",
    "$$\n",
    "p(B_iB_j|I) = \\delta_{ij} = \\begin{cases}1,&\\;\\text{if}\\;i=j,\\\\ 0,&\\;\\text{otherwise}.\\end{cases}\n",
    "$$\n",
    "Then, for any logical sentence $A$ we have:\n",
    "$$\n",
    "p(A|I) = \\sum_{i=1}^n p(AB_i|I) = \\sum_{i=1}^n p(A|B_i,I)p(B_i|I).\n",
    "$$\n",
    "\n",
    "Again, this requires a bit of meditation.\n",
    "You take any logical sentence A and set of exclusive but exhaustive possibilities $B_1,\\dots,B_n$ and you break down the probability of $A$ in terms of the probabilities of the $B_i$'s.\n",
    "The Venn diagrams helps to understand the situation:\n",
    "![Venn sum rule](venn_sum_rule.png)\n",
    "\n",
    "The sum rule can be trivially proved by induction using only the obvious rule and the product rule.\n",
    "It is instructive to go through the proof.\n",
    "For $n=2$ we have:\n",
    "$$\n",
    "\\begin{split}\n",
    "p(A|I) &=& p(A\\;\\text{and}\\;(B_1\\;\\text{or}\\;B_2)|I)\\\\\n",
    "&=& p\\left((A\\;\\text{and}\\;B_1)\\;\\text{or}\\;(A\\;\\text{and}\\;B_2)|I\\right)\\\\\n",
    "&=& p(A\\;\\text{and}\\;B_1|I) + p(A\\;\\text{and}\\;B_2|I) - p\\left((A\\;\\text{and}\\;B_1)\\;\\text{or}\\;(A\\;\\text{and}\\;B_2)|I\\right)\\\\\n",
    "&=& p(AB_1|I) + p(AB_2|I) - p(AB_1B_2|I)\\\\\n",
    "&=& p(AB_1|I) + p(AB_2|I),\n",
    "\\end{split}\n",
    "$$\n",
    "because \n",
    "$$\n",
    "p(AB_1B_2|I) = p(B_1B_2|I)p(A|I) \\le p(B_1B_2|I) = 0.\n",
    "$$\n",
    "And then, assume that it holds for $n$, you can easily show that it also holds for $n+1$ completing the proof."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Drawing balls from a box without replacement (2/3)\n",
    "\n",
    "Let us consider the probability of getting a red ball in the second draw without observing in the first draw $p(B_1|I)$.\n",
    "We have two possibilities for the first draw.\n",
    "We either got a blue ball (B_1 is true) or we got a red ball (R_1 is true).\n",
    "In other words $B_1$ and $R_1$ cover all possibilities and are mutually exclusive.\n",
    "We can use the sum rule:\n",
    "$$\n",
    "\\begin{split}\n",
    "p(R_2|I) &=& p(R_2|B_1,I)p(B_1|I) + p(R_2|R_1,I)p(R_1|I)\\\\\n",
    "&=& \\frac{2}{3}\\frac{2}{5} + \\frac{5}{9}\\frac{3}{5}\\\\\n",
    "&=& 0.6.\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Let's verify this using a Monte Carlo experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_balls_no_replacement(k, num_red, num_blue):\n",
    "    \"\"\"\n",
    "    Draw k balls from an urn without replacement.\n",
    "    \n",
    "    :param k: the number of balls to draw\n",
    "    :param num_red: the number of red balls\n",
    "    :param num_blue: the number of blue balls\n",
    "    :return: a k-element list of 'R' and 'B' depending on the draw.\n",
    "            If k is greater than num_red + num_blue, all the balls are returned.\n",
    "    \"\"\"\n",
    "    balls = ['R'] * num_red + ['B'] * num_blue\n",
    "    np.random.shuffle(balls)\n",
    "    drawn_balls = []\n",
    "    for i in range(k):\n",
    "        # Pick a ball at random\n",
    "        ball_id = np.random.randint(0, len(balls))\n",
    "        ball = balls[ball_id]\n",
    "        # Remove it from the urn\n",
    "        balls = balls[:ball_id] + balls[ball_id+1:]\n",
    "        # Put it in the drawn balls\n",
    "        drawn_balls.append(ball)\n",
    "    return drawn_balls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['B', 'R', 'R']"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# You can use it like this:\n",
    "draw_balls_no_replacement(3, 6, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of red on second draw: 0.60\n"
     ]
    }
   ],
   "source": [
    "# Now let's draw two balls a few thousand times and count how many times we get a red one\n",
    "# in the **second** draw\n",
    "num_samples = 100000\n",
    "running_sum_of_R2 = 0\n",
    "for n in range(num_samples):\n",
    "    balls = draw_balls_no_replacement(2, 6, 4)\n",
    "    if balls[1] == 'R':\n",
    "        running_sum_of_R2 += 1\n",
    "p_R2 = running_sum_of_R2 / num_samples\n",
    "print('Probability of red on second draw: {0:1.2f}'.format(p_R2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which matches exactly the theoretical result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ What is the probability of a blue ball on the second draw $p(B_2|I)$.\n",
    "Find the theoretical result and then modify the Monte Carlo above to verify it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Drawing balls from a box without replacement (3/3)\n",
    "\n",
    "If you paid close attention, in all our examples the conditioning we did followed the causal links.\n",
    "For instance, in the urn example we where writing $p(R_2|B_1,I)$ for the probability of getting a red ball in the second draw after having observed the blue ball in the first draw.\n",
    "This is the uncertainty propagation problem.\n",
    "However, conditioning on stuff **does not have to follow the causal links**.\n",
    "It is completely legitimate to ask what is the probability of a blue ball in the first draw given that you have observed that the result of the second draw is a red ball.\n",
    "The situation is visualized in the following graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Urn4 Pages: 1 -->\n",
       "<svg width=\"227pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 226.82 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\n",
       "<title>Urn4</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-184 222.82,-184 222.82,4 -4,4\"/>\n",
       "<!-- reds -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>reds</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"46.94\" cy=\"-162\" rx=\"46.89\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"46.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\"># red balls</text>\n",
       "</g>\n",
       "<!-- first -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>first</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"104.94\" cy=\"-90\" rx=\"40.13\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">1st draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;first -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>reds&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M60.69,-144.41C67.91,-135.69 76.9,-124.85 84.89,-115.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"87.61,-117.4 91.3,-107.47 82.22,-112.94 87.61,-117.4\"/>\n",
       "</g>\n",
       "<!-- second -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>second</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"104.94\" cy=\"-18\" rx=\"43.02\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.94\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">2nd draw</text>\n",
       "</g>\n",
       "<!-- reds&#45;&gt;second -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>reds&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.54,-143.77C44.73,-125.21 45.48,-95.22 55.94,-72 61.04,-60.69 69.53,-50.2 78.01,-41.59\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"80.58,-43.98 85.39,-34.53 75.74,-38.92 80.58,-43.98\"/>\n",
       "</g>\n",
       "<!-- blues -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>blues</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"167.94\" cy=\"-162\" rx=\"50.76\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"167.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\"># blue balls</text>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;first -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>blues&#45;&gt;first</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M153.01,-144.41C144.98,-135.48 134.92,-124.31 126.09,-114.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"128.63,-112.08 119.33,-106.99 123.42,-116.76 128.63,-112.08\"/>\n",
       "</g>\n",
       "<!-- blues&#45;&gt;second -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>blues&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M168.15,-143.72C167.73,-125.12 165.13,-95.09 153.94,-72 148.45,-60.66 139.71,-50.05 131.13,-41.35\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"133.35,-38.63 123.71,-34.22 128.5,-43.68 133.35,-38.63\"/>\n",
       "</g>\n",
       "<!-- first&#45;&gt;second -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>first&#45;&gt;second</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M104.94,-71.7C104.94,-63.98 104.94,-54.71 104.94,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"108.44,-46.1 104.94,-36.1 101.44,-46.1 108.44,-46.1\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a22fc4890>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gu4 = Digraph('Urn4')\n",
    "gu4.node('reds', label='# red balls', style='filled')\n",
    "gu4.node('blues', label='# blue balls', style='filled')\n",
    "gu4.node('first', label='1st draw')\n",
    "gu4.node('second', label='2nd draw', style='filled')\n",
    "gu4.edge('reds', 'first')\n",
    "gu4.edge('blues', 'first')\n",
    "gu4.edge('first', 'second')\n",
    "gu4.edge('reds', 'second')\n",
    "gu4.edge('blues', 'second')\n",
    "gu4.render('urn4_graph', format='png')\n",
    "gu4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, you can write down the mathematical expression $p(B_1|R_2,I)$.\n",
    "This does not mean that $R_2$ is causing $B_1$.\n",
    "What happens here is that observing $R_2$ changes your state of knowledge about $B_1$.\n",
    "This is an example of information flowing in the reverse order of a causal link and a quintessential example of the inverse problem.\n",
    "Let's solve it analytically:\n",
    "$$\n",
    "\\begin{split}\n",
    "p(B_1|R_2,I) &=& \\frac{p(B_1,R_2|I)}{p(R_2|I)}\\\\\n",
    "&=& \\frac{\\frac{4}{15}}{0.6}\\\\\n",
    "&\\approx& 0.44.\n",
    "\\end{split}\n",
    "$$\n",
    "This is greater than the probability of drawing a blue ball in the first place, $p(B_1|I) = 0.4$.\n",
    "Does this make sense?\n",
    "Yes it does!\n",
    "Here is how you should think:\n",
    "+ You draw a ball without seeing it and you put in a box.\n",
    "+ You draw the second ball and you see that it is a red one.\n",
    "+ This means that this particular red ball was not picked in the first draw.\n",
    "+ So, it is as if in the first draw you had one less red to worry about which increases the probability of a blue.\n",
    "+ So, it is as if you had 5 red balls and 4 blue balls giving you a probability of blue $\\frac{4}{9}\\approx 0.44$.\n",
    "\n",
    "This is amazing!\n",
    "It agrees perfectly with the prediction of the product rule.\n",
    "This was one of our desiderata (if you compute something in two different ways you should get the same result).\n",
    "You can rest assured that as soon as you use the product rule and the sum rule and logic, it is impossible to get the wrong answer.\n",
    "That is, if you can actually carry out the computations.\n",
    "\n",
    "Let's also verify the result with a Monte Carlo simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of blue on first draw: 0.447\n"
     ]
    }
   ],
   "source": [
    "num_samples = 100000\n",
    "running_sum_of_R2 = 0\n",
    "running_sum_of_B1_given_R2 = 0\n",
    "for n in range(num_samples):\n",
    "    balls = draw_balls_no_replacement(2, 6, 4)\n",
    "    if balls[1] == 'R':\n",
    "        running_sum_of_R2 += 1\n",
    "        if balls[0] == 'B':\n",
    "            running_sum_of_B1_given_R2 += 1\n",
    "p_B1_given_R2 = running_sum_of_B1_given_R2 / running_sum_of_R2\n",
    "print('Probability of blue on first draw: {0:1.3f}'.format(p_B1_given_R2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Work out theoretically $p(R_1|B_2,I)$.\n",
    "+ Verify your result by modifying the Monte Carlo procedure above."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  },
  "widgets": {
   "state": {
    "968a18c908c6495081a9936faebe0473": {
     "views": [
      {
       "cell_index": 11
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
